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SUMMARY 


Simulation  of  multi-phase,  turbulent  reacting  flow  is  in  itself  a  very  complex  task 
but  when  such  flows  occur  in  the  presence  of  strong,  unsteady  shocks  additional  com¬ 
plexity  can  arise.  Shock  interactions  with  shear  turbulence  can  change  turbulent 
structures  and  shock  induced  heating  can  trigger  ignition,  combustion  and  turbulent 
flame  propagation.  There  are  many  practical  applications  of  interest  to  US  Air  Force, 
such  as  scramjets,  pulse  detonation  engines,  multi-phase  detonation  (explosives  con¬ 
taining  inert  or  reactive  particles),  and  detonation  interaction  with  reactive  mixtures 
where  multi-phase  (gas-liquid,  gas-solid)  physics  couple  with  shock-turbulence-flame 
interactions.  Simulation  of  such  flows  (either  using  direct  or  large-eddy  simulation 
approach)  requires  a  numerical  strategy  that  can  simultaneously  capture  strong  shock 
motion  and  turbulence-flame  interactions  occurring  over  a  wide  range  of  spatial  and 
temporal  scales.  In  this  research,  a  new  and  an  efficient  large-eddy  simulation  (LES) 
strategy  has  been  developed  to  investigate  turbulent  flows  in  a  high-speed,  compress¬ 
ible  environment.  A  new  numerical  algorithm  has  been  developed  that  permits  a 
proper  capture  of  the  discontinuities  that  arise  at  supersonic  speeds  using  a  hybrid 
Riemann  solver  that  shows  reduced  sensitivity  to  classical  solver  instabilities.  Since  an 
upwind  shock-capturing  schemes  can  dissipate  turbulent  features  the  current  hybrid 
approach  uses  a  low  dissipation  algorithm  in  regions  of  high  shear  without  shocks. 
In  order  to  simulate  high  Reynolds  number  flows,  large-eddy  simulation  (LES)  using 
a  new  compressible  subgrid  closure  has  been  developed  and  validated.  Direct  simu¬ 
lations  of  shock/turbulence  interactions  for  Mach  numbers  ranging  from  1.29  to  3.0 
have  been  used  to  identify  the  most  important  terms  in  the  closure  strategy  adopted 
for  shock-containing  flows.  The  models  for  these  terms  and  the  dynamic  evaluation 
of  the  closure  coefficients  as  a  part  of  the  solution  demonstrated  in  this  study  lead  to 
a  very  general  methodology  with  no  ad  hoc  constants.  This  closure  model  has  been 
validated  by  performing  Large  Eddy  Simulations  of  the  shock  /  turbulence  flow  config¬ 
urations  as  well  of  other  flow  configurations,  such  as  Richtmyer-Meshkov  instabilities, 
regular  and  irregular  gaseous  detonations,  and  mixing-augmentation  configurations 
based  on  shock  /  turbulent  shear  interactions. 


CHAPTER  I 


INTRODUCTION 


A  class  of  problems  of  both  great  fundamental  interest  and  practical  relevance  is  in 
the  held  of  highly  compressible  turbulent  hows  of  multi-huids.  Typical  challenging 
problems  are  in  multi-phase  detonations  (these  can  involve  gaseous  detonation  wave 
interacting  with  two-phase  reactive  mixture  or  detonation  products  containing  re¬ 
active  or  inert  particles),  and  in  strong  shock  wave  propagation  in  turbulent  media 
followed  by  shock  induced  ignition  and  combustion.  Shock  interaction  with  turbu¬ 
lence  and/or  hames  have  many  practical  applications  for  US  Air  Force  in  applications 
such  as  scramjets,  pulse-detonation  engines  (PDF),  stage  separation,  supersonic  cav¬ 
ity  oscillations,  hypersonic  aerodynamics,  detonation  induced  structural  destruction, 
detonation  induced  destruction  of  chemical  and  biological  agents,  etc. 

In  shock-turbulence  interactions,  post-shock  and/or  post-detonation  turbulence 
evolution  is  of  particular  interest  from  a  fundamental  perspective  as  well  since  some 
of  the  turbulent  physics  may  not  follow  conventional  well  known  scaling  laws.  For 
example,  extensive  studies  of  Richtmyer-Meshkov  instability  (RMI)  are  currently  un¬ 
derway,  both  experimentally  [55,  2,  9]  and  numerically  [11,  28,  49]  since  RMI  is 
a  fundamental  instability  that  occurs  when  a  passing  shock  wave  interacts  with  a 
perturbed  density  interface  separating  gases  of  diherent  densities  [5,  57].  Vorticity 
generation  occurs  around  the  interface  and  during  subsequent  evolution  develops  into 
many  scales  of  turbulent  like  structures.  If  this  turbulent  held  is  re-shocked,  addi¬ 
tional  scales  of  turbulence  have  been  shown  to  be  generated.  RMI  and  reshocked-RMI 
[59,  60,  49]  are  very  relevant  for  core-collapse  supernova  dynamics,  but  can  also  occur 
in  the  detonation  scenarios  involving  reacting  hows,  including  supersonic  combustion 
ramjets  (SCRAMJETS)  and  pulse  detonation  engines.  Shock  induced  mixing  changes 
associated  with  RMI  and  re-shocked  RMI  are  therefore,  of  considerable  interest  from 
practical  point  of  view.  There  are  also  other  fundamental  unresolved  issues  regard¬ 
ing  RMI.  For  example,  recent  studies  in  re-shocked  RMI  suggest  that  the  scales  of 
turbulence  in  the  post  shock  region  may  not  strictly  follow  the  classical  Kolmogorov 
—5/3  inertial  range  scaling  [59,  44].  If  so,  then  this  type  of  turbulence  may  require 
diherent  scaling  arguments,  and  hence  new  types  of  subgrid  closure(s). 

There  is  very  little  experimental  data  (even  when  available  the  data  is  sparse  and 
not  time-resolved)  and  these  multi-physics  problems  are  also  inherently  difficult  to 
solve  due  to  the  very  large  range  of  temporal  and  spatial  scales  involved.  Numerical 
approach  must  capture  strong  moving  shocks  and  hne-scales  of  turbulence.  Con¬ 
ventional  shock  capturing  schemes  are  too  dissipative  for  this  purpose.  As  a  result, 
new  high-order  schemes,  such  as  the  ninth-order  weighted  essentially  non-oscillatory 
(WENO)  shock-capturing  schemes  have  been  developed  for  direct  numerical  simula¬ 
tions  (DNS).  Most  of  the  past  RMI  simulations  (primarily  using  DNS)  using  such 
schemes  have  been  restricted  to  2D  [23,  49,  33]  and  only  a  few  studies  have  reported 
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using  full  3D  DNS  or  LES  approaches  [11,  28].  Furthermore,  all  these  studies  are 
in  very  simplihed  (planar)  shock-tube  type  geometries.  However,  for  practical  appli¬ 
cation  to  SCRAMJET  or  PDE  type  applications,  the  geometrical  complexities  and 
test  conditions  are  such  that  DNS  and  even  very  high  order  schemes  are  not  com¬ 
putationally  practical.  On  the  other  hand,  large-eddy  simulation  (LES)  may  be  a 
viable  alternative  since  it  has  already  shown  potential  for  subsonic  problems  [31,  41]. 
However,  in  order  to  conduct  LES  of  supersonic  and  shock-dominated  flows  new  al¬ 
gorithms  and  subgrid  closures  for  highly  compressible  flow  have  to  be  developed  and 
validated.  LES  studies  of  compressible  mixing,  shock-shear  interactions,  detonations 
and  RMI  have  just  started  to  appear  [28,  22,  25,  24]  but  there  are  still  many  issues  to 
be  addressed.  For  example,  in  addition  to  the  scales  associated  with  the  shock  struc¬ 
ture  and  hne-scale  turbulence,  particle  (e.g.,  liquid  droplets)  motion,  fuel-air  mixing 
and  hnite-rate  reactions  have  to  be  included  for  some  applications.  Chemical  kinetics 
are  usually  very  stiff  requiring  very  small  time-steps  to  properly  resolve  turbulence- 
chemistry-shock  interactions.  Implicit  schemes  have  only  limited  functionality  (since 
there  is  a  limit  beyond  which  the  time  scales  of  interactions  have  to  be  resolved),  and 
also  have  problems  scaling  in  massively  parallel  systems. 

Thus,  there  are  still  many  challenges  to  study  of  shock-turbulence-flame  (STF) 
interaction  in  both  canonical  and  practical  systems.  Some  of  these  challenges  have 
been  addressed  over  the  last  few  years  at  Georgia  Tech,  and  has  resulted  in  a  new 
multi-scale,  multi-phase  solver  that  can  simulate  STF  interactions  with  good  accuracy. 
This  capability  was  developed  in  part  using  support  from  AFOSR  and  the  code  has 
been  delivered  to  Eglin  AFB  to  support  their  in-house  (classihed)  applications.  In 
fact,  the  applications  of  this  code  in  Eglin  AFB  to  problems  of  AF  practical  interest 
demonstrates  our  ability  to  transition  from  fundamental  studies  conducted  in  GT  to 
more  complex  practical  problems  of  AF  interest. 

This  report  summarized  the  key  results  obtained  in  the  last  hscal  year  of  this 
research  effort.  Since  this  is  also  its  hnal  year,  some  key  issues  and  developments 
reported  in  previous  reports  have  been  included  for  completeness.  Further  results 
and  analysis  are  reported  in  cited  publications. 
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CHAPTER  II 


SUMMARY  OF  RESULTS 


To  establish  a  validated  simulation  capability  for  STF  interaction  studies,  we  leverage 
a  code  called  LESLIE3D  (Large-Eddy  Simulation  using  Linear  Eddy  in  3D),  origi¬ 
nally  developed  to  study  turbulent  two-phase  combustion  in  subsonic  flows  as  in  gas 
turbine  and  ramjet  combustors.  This  code  is  a  multi-block,  structured  grid  solver  of 
the  fully  compressible  Navier-Stokes  equations  using  a  hnite-volume  scheme  that  is 
second-order  accurate  in  time  and  up  to  fourth-order  accurate  in  space.  In  addition 
to  LES  the  code  can  be  used  to  perform  DNS  and  has  been  extensively  validated  in 
many  non- reacting  and  reacting  flow  apphcations(see  [41]  and  other  cited  references). 

To  extend  this  code  to  handle  shock  and  detonation  physics  two  key  challenges  have 
been  continuously  addressed  and  successfully  resolved  in  the  current  effort.  The  hrst 
challenge  consists  in  developing  a  numerical  approach  that  adapts  locally  to  the  un¬ 
steady  shock  motion  and/or  compressive  regions  without  contaminating  or  dissipating 
the  turbulent  and/or  flame  structures.  This  requires  developing  and  validating  a  new 
hybrid  approach  within  the  generalized  multi-block  approach  that  can  achieve  simul¬ 
taneously,  very  low  dissipation  in  smooth  flow  regions  and  proper  shock  capturing  in 
the  regions  of  high  compressibility.  The  second  challenge  is  to  develop  improved  LES 
sub-grid  closure  model  to  include  the  effects  of  high  compressibility  on  the  turbulence 
evolution.  A  new  localized  dynamic  sub-grid  closure  that  can  handle  some  of  the 
underlying  physics  of  STF  interactions  has  been  developed  and  validated. 

In  the  following,  we  briefly  summarize  the  key  progress  and  success  of  the  effort 
conducted  over  the  last  3  years.  Some  of  these  results  and  application  of  the  solver 
have  been  reported  in  publications  and  more  detailed  papers  are  currently  being 
written.  Here,  we  list  the  publications  past  and  planned  for  completeness: 

Papers  directly  related  to  research  under  the  current  AFOSR  project: 

1.  Fryxell,  B.  and  Menon,  S.,  “Hybrid  Simulations  of  Richtmyer-Meshkov  Insta¬ 
bility,”  AIAA-2005-0314,  43rd  AIAA  Aerospace  Sciences  Meeting,  Reno,  NV, 
January,  2005. 

2.  Genin,  F.,  Fryxell,  B.,  and  Menon,  S.,  “Simulation  of  Detonation  Propagation  in 
Turbulent  Gas-Solid  Reactive  Mixtures,’”  AIAA-2005-3967,  41st  AIAA/ASME/ASEE 
Joint  Propulsion  Gonference,  Tucson,  AZ,  July  11-14,  2005. 

3.  Genin  F.,  Fryxell,  B.  and  Menon,  S.,  “Hybrid  Large-Eddy  Simulation  of  Det¬ 
onation  in  Reactive  Mixtures,”  20th  International  Gonference  on  Detonations, 
Explosions  and  Shock  Waves,  Montreal  GA,  August  1-4,  2005. 

4.  Menon,  S.,  “A  Multi-Scale  Simulation  Methodology  for  Gompressible  Flows 
and  Turbulent  Gombustion”  Proceedings  of  the  Third  European  Gombustion 
Meeting,  EGM2007,  Ghania,  Grete,  April  11-13,  2007. 
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5.  Fryxell,  B.  and  Menon,  S.,  “Simulation  of  Richtmyer-Meshkov  Instability,”  to 
appear  in  Physica  Scripta,  2008. 

6.  Genin,  F.  and  Menon,  S.,  “Simulation  of  Interaction  of  Shocks  with  Turbulent 
Shear  Layers  Using  Large-Eddy  Simulations,”  Computer  and  Fluids,  2008  (to 
be  submitted). 

7.  Genin,  F.  and  Menon,  S.,  “Simulation  of  Turbulent  Mixing  Behind  a  Strut  Injec¬ 
tor  in  Supersonic  Flow,”  invited  paper,  AIAA-2009-0132,  47th  AIAA  Aerospace 
Sciences  Meeting,  January  2009. 

8.  Ukai,  S.,  Genin,  F.,  Srinivasan,  S.  and  Menon,  S.,  “Large-Eddy  Simulation 
of  Re-shocked  Richtmyer-Meshkov  Instability,”  AIAA-2009-0944,  47th  AIAA 
Aerospace  Sciences  Gonference,  January  2009. 

9.  Gottiparthi,  K.,  Genin,  F.,  Srinivasan,  S.  and  Menon,  S.,  “Simulation  of  Gellular 
Detonation  Structures  in  Ethylene- Air  Mixtures,”  AIAA-2009-0437,  47th  AIAA 
Aerospace  Sciences  Gonference,  January  2009. 

Transition  of  the  numerical  and  subgrid  algorithms  developed  in  this  research  to 
other  research  is  a  key  demonstration  of  the  validity  and  the  relevance  of  the  new 
approaches.  In  this  regard,  we  have  been  very  successful  with  many  publications  (for 
research  funded  by  ONR,  DTRA,  AFOSR  and  NASA)  that  employed  the  algorithms 
developed  under  the  current  AFOSR  project.  Papers  that  uses  the  hybrid  solver 
developed  under  the  AFOSR  project: 

1.  Miki,  K.  and  Menon,  S.,  “Local  Dynamic  Subgrid  Glosure  for  Gompressible 
MUD  Turbulence  Simulation,”  AIAA-2006-2891,  37th  AIAA  Plasma-dynamics 
and  Lasers  Gonference,  San  Francisco,  GA,  5-8  June  2006. 

2.  Schulz,  J.,  Miki,  K.,  and  Menon,  S.,  “Simulation  of  MUD  Turbulence-Ghemistry 
Interaction  in  Supersonic  Flow,”  AIAA  Paper  2006-2894,  36th  AIAA  Plasma- 
dynamics  and  Lasers  Gonf.,  San  Francisco,  5-8  June,  2006. 

3.  Miki,  K.,  Schulz,  J.  and  Menon,  S.,  “Large-Eddy  Simulation  of  Supersonic 
Plasma  Flow  over  a  Backward  Facing  Step,”  Proceedings  of  the  Fifth  Turbulent 
Shear  Flow  Phenomenon  (TSFP-5),  pp.  985-990,  2007. 

4.  Balakrishanan,  K.  and  Menon,  S.,  “Simulation  of  Detonation  Propagation  Through 
Solid  Particle  Glouds  and  Subsequent  Dispersion  of  Particles,”  AIAA-2008-4689, 
44th  AIAA  Joint  Propulsion  Gonference,  Hartford,  GT,  July  21-23,  2008. 

5.  Balakrishanan,  K.,  Genin,  F.,  Menon,  S.,  and  Nance,  D.,  “Numerical  Study 
of  Shock  Overpressure  and  Propagation  Gharacteristics  from  Detonations  in 
Homogeneous  Explosives  Nitromethane,  HMX,  and  TNT,”  2008  (under  prepa¬ 
ration). 
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6.  Balakrishanan,  K.,  Genin,  F.,  Menon,  S.,  and  Nance,  D.,  “Numerical  Study 
of  Shock  Overpressure  and  Particle  Motion  from  Detonations  in  Heterogeneous 
Explosives  containing  Dense  Inert  Particles,”  2008  (under  preparation). 

7.  Miki,  K.,  Schulz,  J.  and  Menon,  S.,  “Large-Eddy  Simulation  of  Non-Equilibrium 
Plasma- Assisted  Combustion  in  Supersonic  Flow,”  AIAA-2008-5072,  44th  AIAA 
Joint  Propulsion  Conference,  Hartford,  CT,  July  21-23,  2008. 

8.  Miki,  K.,  Schulz,  J.  and  Menon,  S.,  “Large-Eddy  Simulation  of  Equilibrium 
Plasma- Assisted  Combustion  in  Supersonic  Flow,”  Proceedings  of  the  Combus¬ 
tion  Institute,  Vol.  32,  2008  (to  appear). 

9.  Miki,  K.,  and  Menon,  S.,  “Localized  Dynamic  Subgrid  Closure  for  Simulation 
of  Magnetohydrodynamics  Turbulence,”  Physics  of  Plasma,  Vol.  15,  072306, 
2008. 

2.1  New  Hybrid  Framework  for  Shock- Turbulence- 
Flame  Interactions 

A  major  problem  for  DNS  and  LES  using  algorithms  developed  for  smooth  flows  is 
that  they  produce  unphysical  oscillations  in  regions  containing  shocks  and  contact 
discontinuities,  and  require  very  rehne  resolution  in  the  discontinuity  region  [34]. 
On  the  other  hand,  shock-capturing  schemes  have  not  shown  potential  for  resolving 
unsteady  shear  flows  in  conjunction  with  shock  waves,  unless  the  resolution  is  high 
enough  so  that  the  dissipation  necessary  to  stabilize  discontinuities  in  typical  shock¬ 
capturing  schemes  do  not  produce  an  unphysical  rate  of  decay  of  turbulent  features. 

In  order  to  capture  strong  unsteady  shock  motion  and  also  the  important  scales  of 
turbulence,  two  key  capabilities  were  developed  and  validated.  For  shock  capturing 
we  initially  (2005-2006)  implemented  Piecewise  Parabolized  Method  (PPM),  which 
is  nominally  0(3)  away  from  shocks  but  is  limited  to  uniform  or  weakly  stretching 
Cartesian  grids.  In  order  to  simulate  shocks  in  more  complex  flows,  we  developed  an 
alternate  high-order  flux  difference  splitting  method  in  2006-2007  within  the  frame¬ 
work  of  our  multi-block  structured  flow  solver.  A  Riemann  problem  is  constructed  at 
the  interface  separating  two  cells,  using  high-order  MUSCL  type  interpolation  (with 
a  Monotonized  Central  limiter  to  prevent  non-monotonic  behavior)  and  solved  us¬ 
ing  an  adaptive  approximate  Riemann  solver.  This  adaptive  approximate  Riemann 
solver  is  based  on  the  Harten-Lax-Leer  (HLL)  family  of  Riemann  solvers  [27].  In 
particular,  we  combine  the  contact  preserving  HLLC  with  a  non-contact  preserving 
HLLE  in  order  to  avoid  parasitic  oscillation  in  the  post  shock  regions.  Detailed  vali¬ 
dation  of  this  combined  HLLC/E  approach  has  been  conducted  and  demonstrated  its 
capability  to  achieve  proper  shock  and  discontinuity  capturing  with  minimal  insta¬ 
bilities.  Switching  between  HLLC  and  HLLE  is  based  on  local  sensors  that  requires 
no  ad  hoc  adjustment  from  the  user.  This  sensor  uses  both  pressure  ratio  and  ve¬ 
locity  divergence.  The  resulting  hybrid  Riemann  solver  differs  from  earlier  hybrid 
solvers  [45,  46]  in  the  sense  that,  instead  of  applying  HLLE  everywhere  in  the  shock 
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thickness,  HLLC/E  uses  HLLE  only  in  the  direction  transverse  to  the  detected  shock. 
This  limits  the  instabilities  in  the  shock  region  and  also  limits  the  amount  of  extra 
dissipation  introduced  by  HLLE.  As  shown  in  our  studies  this  sensor  is  very  robust, 
thanks  to  the  uniqueness  of  the  shock  characteristics  and  has  required  no  ad  hoc  user 
input. 

The  overall  accuracy  of  the  HLLC/E  is  second  order,  although  viscous  terms  are 
still  computed  using  fourth  order  accuracy.  The  HLLC/E  model  is  implemented 
within  our  full  3D  generalized  multi-block,  multi-phase  solvers  LESLIE3D  and  there¬ 
fore,  it  can  exploit  all  the  in-built  capabilities  of  the  code. 

To  simulate  STF  interactions,  another  type  of  hybrid  approach  is  developed, 
whereby,  the  0(4)  smooth  flow  solver  in  LESLIE3D  is  used  in  regions  with  no  shocks, 
and  the  HLLC/E  approach  is  used  in  regions  of  strong  discontinuities.  The  switch 
[22]  employs  detecting  both  pressure  and  density  gradients  and  has  been  successfully 
applied  to  practical  turbulent  problems  [25,  24].  The  regions  of  strong  discontinu¬ 
ities  are  detected  by  a  measure  of  the  smoothness  of  the  solution  locally  so  that  the 
switching  between  the  two  schemes  is  natural,  and  again  requires  no  user  input. 

Thus,  we  now  have  a  baseline  generalized  solver  that  employs  two  layers  of  hybrid, 
locally  adapting  algorithms:  (a)  a  shock  capturing  scheme  switching  to  0(4)  smooth 
flow  solvers  away  from  regions  of  strong  discontinuities,  and  (b)  within  the  regions  of 
strong  discontinuities  a  switch  between  contact  discontinuity  preserving  (HLLC)  and 
non-contact  preserving  (HLLE)  so  that  both  shocks  and  contact  discontinuities  are 
captured  in  full  3D  with  minimal  instabilities.  Application  to  many  canonical  and 
complex  flows  have  demonstrated  the  ability  of  this  new  solver.  Results  for  some  of 
the  validation  studies  are  presented  later. 

2.2  New  Sub-Grid  Closure  for  Compressible  Tur¬ 
bulent  Flows 

The  method  described  above  is  used  to  solve  the  Navier-Stokes  equations  in  general 
but  in  highly  turbulent  flows,  the  grid  resolution  is  insufficient  to  do  direct  numer¬ 
ical  simulation  (DNS)  and  therefore,  we  have  to  resort  to  large-eddy  simulations 
(LES).  LES  capability  in  LESLIE3D  for  subsonic  reacting  flows  is  well  established 
and  is  based  on  a  localized  dynamic  closure  for  the  sub-grid  turbulent  kinetic  energy 
(LDKM).  The  LDKM  approach  allows  the  model  coefficients  to  evolve  as  a  part  of  the 
solution,  and  thus,  there  are  no  model  constants  to  adjust.  This  provides  a  generality 
critical  to  study  complex  flows. 

To  extend  to  highly  compressible  flows,  additional  sub-grid  terms  (that  are  typ¬ 
ically  ignored  in  low-speed  flows)  have  to  be  re-assessed.  For  example,  the  sub-grid 
enthalpy  flux  is  typically  modeled  using  the  sub-grid  eddy  viscosity  and  a  constant 
turbulent  Prandtl  number  in  low  speed  flows.  However,  in  compressible  flow,  the  con¬ 
stant  Prandtl  number  assumption  is  not  accurate  and  therefore,  a  new  scale-similarity 
dynamic  closure  for  the  turbulent  Prandtl  number  has  been  developed.  In  addition, 
the  diffusion  of  due  to  pressure  fluctuations,  which  plays  an  important  role  during 
the  interaction  of  a  turbulent  held  with  a  shock,  has  the  same  structural  form  as  the 
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diffusion  of  energy  due  to  the  subgrid  energy  fluctuations.  In  a  perfect  gas  these  two 
unclosed  terms  are  proportional,  and  therefore,  a  model  for  this  term  is  also  achieved 
through  the  dynamic  computation  of  the  turbulent  Prandtl  number. 

Finally,  the  evolution  of  the  turbulent  kinetic  energy  in  compressible  flows  is 
dependent  on  the  pressure-dilatation  correlation.  This  physical  phenomenon  transfers 
energy  between  internal  and  kinetic  energy,  and  can  be  either  a  source  or  a  sink  of 
kinetic  energy.  An  earlier  analytical  study  [47]  showed  that  the  pressure-dilatation 
correlation  becomes  important  in  non-equilibrium  flows,  that  is,  when  production  and 
dissipation  of  turbulent  kinetic  energy  are  signihcantly  unbalanced.  A  dependence  on 
the  square  of  the  turbulent  Mach  number  is  also  derived  from  this  work  and  used  in 
the  current  closure. 

In  summary,  the  sub-grid  closure  for  compressible  flow  now  employs  additional 
dynamic  closure  for  the  turbulent  Prandtl  number,  the  kinetic  energy  diffusion  due 
to  pressure  fluctuations  and  the  pressure-dilatation  correlation. 

The  developments  described  here  have  also  permitted  high-hdelity  simulations  of 
detonation  processes.  This  combustion  regime  consists  in  a  high-speed  shock  traveling 
in  a  low-temperature  reactive  mixture.  The  gas  is  compressed  and  heated  through  the 
compression,  and  auto-ignition  is  achieved.  The  heat  released  by  the  reaction  feeds 
back  the  main  wave.  A  good  capture  of  the  phenomenon  requires  a  proper  capture 
of  the  running  shock,  an  accnrate  treatment  of  the  chemical  sonrce  terms,  and  a 
low  dissipation  in  the  post-shock  region,  in  order  to  capture  the  turbulent  structures 
emanating  from  the  triple  points  on  the  detonation  front.  These  capabilities  have 
been  developed  and  validated  in  onr  solver. 

In  the  following  chapters  we  discuss  the  key  results  of  this  study. 
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CHAPTER  III 


FORMULATION  AND  ALGORITHMS 


3.1  Governing  equations 

3.1.1  Filtered  Navier- Stokes  Equations 

The  direct  resolution  of  the  Navier-Stokes  equations  requires  an  increasingly  high 
computational  expense  as  the  flow  Reynolds  number  is  increased.  Except  for  a  few 
canonical  flows,  this  method  cannot  be  adopted  for  the  simulation  of  practical  flows, 
and  so  LES  with  appropriate  subgrid  closure  is  a  viable  option.  The  compressible 
LES  equations  are  obtained  by  spatially  hltering  the  Navier-Stokes  equations  using 
a  top-hat  hlter  (see,  e.g.  [56]),  appropriate  for  hnite  volume  schemes.  If  we  dehne 
/  as  the  spatially  hltered  variable  /,  and  /  as  the  Favre-hltered  variable  /,  dehned 
by  /  =  pf /p,  where  p  represents  the  local  fluid  density,  then  the  LES  equations  for 
mass,  momentum,  total  energy  and  species  conservation  can  be  written  as: 


dt  dxi 

^  ^  -  V  +  1  =  0 


^  ^  +  p)  “i  +  9i  -  “Vy  +  H‘“  +  (t‘‘ 


(3.1) 

(3.2) 

(3.3) 


^  ^  (us, :  + Ur«)  +  y.:r  +  "S' ]  =  0  k  =  l,...,N.  (3.4) 

Here,  p  is  the  density,  (ui)i=i,2,3  is  the  velocity  vector  in  Cartesian  coordinates,  P 
is  the  pressure,  T  is  the  temperature,  and  is  the  k  —  th  species  mass  fraction.  Also, 
Ng  is  the  total  number  of  species.  The  total  energy  of  the  system  is  the  sum  of  the 
internal  energy  e  and  the  kinetic  energy.  As  a  consequence,  the  hltered  total  energy 
is  the  sum  of  the  hltered  internal  energy  e,  the  resolved  kinetic  energy,  {l/2)[uiUj\ 
and  the  resolved  subgrid  kinetic  energy  given  by  k^^"^  =  {l/2)[uiUi  —  UiUi].  A 

perfect  gas  equation  of  state  is  assumed  here  and  neglecting  the  correlation  between 
subgrid  species  and  the  subgrid  temperature  huctuations,  the  equation  of  state  reads: 


(3.5) 

(3.6) 


where  the  mixture  gas  constant  R  is  obtained  using  the  specie  molecular  weights 
{MWk)k=i,...,NS  and  Ru  the  universal  gas  constant,  as: 


R 


k=l,...,NS 


Yk 


Rn 


MWk 


(3.7) 


The  hltered  static  enthalpy  per  unit  mass  h  is  dehned  as  ph  =  pe  +  P  and  for  perfect 
gas  can  be  expressed  as: 

h=  Y.  (3.8) 

k=l,...,NS 

where  Cp^k  =  c^,fc  +  Rk-  The  specihc  heat  at  constant  pressure  k  —  th  species  k,  Cp^k 
are  all  obtained  from  classical  temperature  curve-hts  [26]  as  polynomials  of  order  4: 

Cp,kiT)  =  Y  (2-^) 

i=0,...,4 


h. 


-iT)=  r  c,,^(T)dT^  Y.  — 

•^0  i=l,...,5  * 


(3.10) 


The  hltered  viscous  stress  tensor,  Tj^,  and  the  hltered  heat-hux  vector,  g*  are  approx¬ 
imated  as: 


Yj  2p.(T')  ^Sij  ^Skk^ij 


Qi  =  -k{T)—  +  p  Yj  hkXkVi^k 

*  fc=l,...,V5 

In  the  previous  expression,  Sij  represents  the  rate  of  strain,  given  as: 


(3.11) 

(3.12) 


’  -  1  ^  I 

2  dxi ) 


(3.13) 


Empirical  temperature-dependent  laws  are  used  for  the  evaluation  of  the  viscosity. 
Following  [35],  the  viscosity  is  computed  using  a  power-law  temperature  relation,  with 
an  exponent  of  0.76  for  the  case  of  shock  /  isotropic  tnrbulence  interaction,  and  a 
Sutherland  law  is  used  for  shock  /  shear  interaction  problem.  In  all  cases,  the  thermal 
condnctivity  is  obtained  by  assnming  a  constant  Prandtl  nnmber  (P^  =  Cpp,/ n  =  0.72). 
Finally,  Id,fc  is  the  hltered  dihnsion  velocity  for  k  —  th  species  and  is  modeled  using 
Pick’s  law: 


i,k 


-Dk  dYk 
Yk  dxi 


(3.14) 


where  Dk  is  the  D^-species  dihnsion  coefficient. 

All  the  subgrid-scale  terms,  denoted  with  a  sgs  superscript  require  specihc  mod¬ 
eling  to  close  the  above  system  of  LES  equations.  These  terms  are  dehned  as: 


(3.15) 
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=  P  (yEui  -  Eu^  +  [uiP  -  UiP)  (3.16) 

-  UjfPj)  (3.17) 

F7f  =  (3.18) 

The  closure  strategy  to  model  the  terms  denoted  as  subgrid-scale  will  be  presented 

next. 


3.1.2  Subgrid  Closure  Model 


An  eddy-viscosity  type  closure  is  adopted  in  this  study.  The  eddy  viscosity,  Vf,  is 
evaluated  using  a  characteristic  length-scale,  provided  by  the  local  grid  size  (A),  and 
a  characteristic  subgrid  velocity,  which  is  obtained  from  the  subgrid  kinetic  energy 
so  that  Vt  =  c^Ay/k^.  The  unclosed  terms  in  the  momentum  equation,  the 
subgrid  stresses  are  then  closed  as: 


sgs  _ 

ij 


‘^P^t  (  Sij  ^ 


(3.19) 


The  two  unclosed  terms  in  the  energy  equation,  and  are  modeled  together 
following: 


(3.20) 


The  subgrid  diffusion  of  species  mass  fractions,  is  also  modeled  using  an  eddy- 
diffusivity  assumption,  as: 


P^t  dYk 
Sct  dxi 


(3.21) 


Thus,  to  close  the  hltered  equations  of  mass,  momentum  and  energy  conservation 
using  the  above  relations  the  local  value  of  must  be  determined.  In  the  present 
study,  a  transport  model  for  the  subgrid  kinetic  energy  is  used.  The  exact  governing 
equation  for  reads: 


Wt 


V  jy.sgs  +  Pk  sgs  p  Pf^sgs  —  Df^sgs 


(3.22) 


where  V^ags  represents  the  diffusion  of  /c®^®,  Vk^gs  is  the  pressure  dilatation  correlation, 
and  Pj^ags  and  D^aga  are  the  production  and  dissipation  of  A;®®®,  respectively.  Their 
exact  expressions  are  given  by: 

Vkaga  =  (^{pK Ui  -pKui  -  UjT-f)  -1-  {uiP  -  UiP)  -  {ufTj  -  Uj  t^) j  (3.23) 


^  ^dui  -^dui 

Vkaga  =  P  —  P  — 

OXi  OXi 

dui 


Pk.aga  —  — T.„“ 


sgs  ^ 


dxi 


(3.24) 

(3.25) 
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(3.26) 


n  _  f  dui  _dui\ 

The  pressure-dilatation  correlation  and  the  diffusion  and  dissipation  of  cannot  be 
already  obtained,  and  require  additional  modeling.  By  analogy  with  the  Kolmogorov 
concept  of  energy  cascade,  the  dissipation  of  subgrid  kinetic  energy  is  evaluated  based 
on  the  characteristic  length-  and  velocity-  scales  as  D/^sgs  =  The 

diffusion  due  to  subgrid  fluctuations  in  kinetic  energy,  subgrid  fluctuations  in  viscous 
stress,  and  subgrid  fluctuations  in  pressure  all  contribute  to  the  global  diffusion  of 
and  each  require  proper  modeling  each.  The  hrst  contribution  (often  referred 
to  as  the  triple  velocity  correlation)  and  the  second  are  modeled  together  using  a 
ksgs  g2;adient  diffusion  approach.  The  diffusion  due  to  pressure  fluctuations  is  often 
neglected  in  incompressible,  or  low-compressibility.  It  has,  however,  been  found  from 
direct  simulations  of  shock  /  turbulence  interaction  (presented  in  last  year’s  report) 
that  it  plays  an  important  role  in  the  /c^^^-evolution  in  shock-containing  flows.  This 
term  isnow  closed  using  an  eddy-viscosity  assumption  as: 

UiP  -  UiP  =  pRiuiT  -  Uif)  (3.27) 

ap  oxi 

So  that  the  global  model  for  the  diffusion  of  subgrid  kinetic  energy  reads: 


V 


ksgs 


A 

dxi 


putR  dT 

dxi  (jp  dxi 


(3.28) 


The  pressure-dilatation  correlation  has  been  the  focus  of  many  studies  in  the 
past,  mostly  from  DNS  or  fundamental  analyses  [1,  58,  48,  47].  Many  of  these  models 
were  unfortunately  obtained  from  DNS  results  with  non  stable  initializations.  The 
analytical  formulation  of  Ristorcelli  [47]  is  performed  in  the  acoustic  limit,  and  is 
adequate  in  the  limit  of  small  turbulent  Mach  numbers.  From  this  analysis,  it  appears 
that  the  pressure-dilatation  correlation  is  found  to  be  an  out-of- equilibrium  turbulent 
phenomenon,  i.e.,  when  the  turbulent  production  is  not  balanced  by  the  turbulent 
dissipation.  The  expression  obtained  by  [47]  reads: 

— -  j  {Pk^gs  —  Dk»9s)  (3.29) 

Uksgs  j 

where  is  the  turbulent  Mach  number  based  on  the  local  value  of  and  the 

local  speed  of  sound.  This  term  becomes  important  when  the  flow  is  both  out  of 
equilibrium  and  the  gradient  Mach  number  [MtSk'^^^ /Dk^gs)  is  large.  We  employ  this 
model  in  our  closure  as  well. 

3.1.3  Localized  dynamic  closure  for  subgrid  model  coefficients 

The  closure  terms  in  the  LES  and  the  equation  for  have  six  closure  coefficients: 
{cu,Cf:,ak,  Prt,apd,  Sct).  Nominal  values  for  these  coefficients  have  been  proposed 


11 


in  the  past  (see,  e.g.  [7,  8]).  In  the  present  study,  a  dynamic  evaluation  of  the 
coefficients  is  performed,  based  on  the  framework  of  the  Localized  Dynamic 
Model  (LDKM)  [32,  30].  This  model  uses  the  experimentally  observed  [38]  similarity 
between  subgrid  stresses  and  Leonard  stresses  to  evaluate  the  closure  coefficient  from 
directly  computable  quantities.  The  similarity  hypothesis  implies  that  the  Leonard 
stress  is  similar  in  nature  to  the  subgrid  stress.  In  particular,  the  trace  of  these 
tensors,  (p)  and  respectively,  are  similar. 

Noting  (/)  the  explicit  application  of  a  top-hat  Liter  of  size  A  (A  >  A)  performed 
on  variable  /,  the  Leonard  stress  can  be  expressed  as: 


{pUj)  {pUj)  \ 

(p)  (p) ) 


(3.30) 


In  the  current  closure,  the  Leonard  stress  at  the  test  filter  level  (which  is  known  explic¬ 
itly  is  related  to  the  subgrid  model  at  the  test  hler  level  by  the  similarity  assumption. 
Thus,  using  the  test-scale  kinetic  energy  (p*®®*  =  l/2£jj)  as: 


C 


-2  (p)  Ci^v^F^A 


1 

3 


I  (P>  k‘“%, 


(3.31) 


so  that  this  relation  has  c,y  as  the  only  remaining  unknown.  Cij  being  a  symmetric 
tensor,  equation  3.31  leaves  us  with  five  independent  relations  for  the  determination 
of  one  unknown.  The  over-determination  is  circumvented  by  using  a  least-square 
method([36])  so  that  Ci,  is  computed  as: 


2MijMij 


where: 


M^j  =  Vk^  (p)  A 


1  ^  ^ 


The  exact  governing  equation  for  P*®®*  can  be  written  as: 


d  d 

-p— p  k  -|~  7: -  (ipHi)  k^  )  =  'Dj^test  Pj^test  —  D j^test 

at  axi  ^  ^ 


(3.32) 


(3.33) 

(3.34) 


(3.35) 


where  the  expression  for  diffusion,  pressure  correlation,  production  and  dissipation  at 
the  test-scale  level  are  fully  expressed  as  functions  of  the  resolved  variables  and  the 
subgrid  stresses.  In  particular,  the  dissipation  is  given  by: 


Potest  —  / 


sgs 


du^ 


dUn 


(3.36) 
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Using  the  similarity  assumption,  the  dissipation  of  can  be  modeled  similar  to 
the  model  for  the  dissipation  of  Detest  =  where  again,  only  the 

closure  coefficient  remains  unknown.  Assuming  that  the  viscosities  can  be  extracted 
from  the  hltering  operation,  the  last  coefficient,  c^,  are  expressed  as: 


Ce 


A 

3" 

(p}(fctest)2 


2 

3 


(3.37) 


where  T^j  is  twice  the  traceless  strain  tensor: 


'^ij  2  I  Sij 


(3.38) 


This  expression  models  both  the  solenoidal  and  dilatational  dissipation  together,  using 
the  same  closure  expression.  This  approximation  remains  acceptable  as  long  as  the 
turbulent  Mach  number  remains  small,  so  that  dilatational  dissipation  remains  very 
small.  It  should  be  noted  that  in  most  applications  of  supersonic  to  hypersonic  flows, 
the  turbulent  Mach  number  does  remain  small,  i.e.,  far  from  unity  [17]. 

In  the  current  study,  the  extension  of  the  LDKM  model  to  compressible  flows 
requires  computing  the  local  turbulent  Prandtl  number  dynamically  using  the  sim¬ 
ilarity  hypothesis.  At  the  test-scale  level,  the  expression  for  the  temperature  and 
velocity  correlation  can  be  expressed  exactly  and  the  over-specified  system  is  again 
solved  using  a  least-square  method,  leading  to  a  final  expression: 


where 


1  /  D  diTli 

^/Pn  =  ,  , 

di  di 

(3.39) 

(3.40) 

„  {kP)  {-pn.)  {t>^) 

(P)  (P)  (P) 

(3.41) 

Similarly,  the  closure  coefficient  for  the  pressure-dilatation  correlation  is  evaluated 
from  the  test-scale  pressure-dilatation  correlation  as: 


^pd 


ac^ 


S)  A^{Pktest 


(3.42) 


where  c  is  the  speed  of  sound. 

The  other  coefficients  Ck  and  Sct  could  also  be  evaluated  dynamically,  if  needed 
using  a  similar  strategy.  However,  the  impact  of  the  ak  is  small  compared  to  other 
terms  in  the  governing  equation  for  /c®^®,  and  is  therefore,  assumed  to  be  constant 
{ak  =  1).  The  turbulent  Schmidt  number  is  taken  to  be  equal  to  0.9.  This  latter 
approximation  is  acceptable  for  the  simulation  of  cases  where  species  play  a  passive 
role. 
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3.1.4  Realizability  conditions 


Vreman  et  al  [56]  showed  that  if  a  positive  semi-dehnite  hlter,  such  as  the  top-hat 
hlter,  is  used  within  an  LES  formulation,  the  subgrid  stresses  have  to  be  positive  semi- 
dehnite.  These  conditions,  referred  to  as  the  realizability  constraints,  were  found  by 
Nelson  and  Menon  [42]  to  be  satished  by  the  LDKM  formulation  for  the  computation 
of  well-resolved  and  fully  turbulent  simulations.  The  strong  and  very  localized  vari¬ 
ations  induced  by  shocks  can  makes  it  difficult  to  satisfy  this  property  over  regions 
of  the  how.  Genin  et  al  [24]  showed  that  enforcing  the  weak  realizability  constraints 
was  benehcial  to  the  LES  modeling  of  shock  /  isotropic  turbulence  interaction  at  low 
Mach  numbers.  The  weak  realizability  constraint  consists  in  enforcing: 


\sgs 
0  12 

2 

< 

sgs  sgs 

Rl  ^2  2 

\sgs 
\  '23 

2 

< 

sgs  sgs 

''’22  ''’33 

(3.43) 

\sgs 
0  13 

2 

< 

sgs  sgs 

Ri  ''’33 

With  the  closure  adopted  in  the  present  study,  these  constraints  are  manipulated  to 
show  that  the  following  relation  should  always  hold: 


Cu  <  ^ 

^  A 


(3.44) 


where 

^  (4-  -  (4-  -  (3-45) 

The  validity  of  the  current  closure  approach  to  these  constraints  will  be  closely  ex¬ 
amined  in  the  context  of  shock  /  turbulence  interaction  for  strong  shocks. 


3.2  Hybrid  Algorithm  for  Shock-Shear  Simulation 

An  accurate  simulation  of  turbulence  in  supersonic  flows  is  made  difficult  by  the  very 
different  scheme  properties  required  within  the  same  simulation.  The  resolution  of 
turbulent  structures  requires  a  smooth-flow  solver  with  limited  numerical  dissipation, 
whereas  the  discontinuities  that  can  arise  from  supersonic  flows  require  a  signihcant 
amount  of  dissipation  in  order  to  be  correctly  captured.  Shock-capturing  method¬ 
ologies  have  been  developed  to  enable  the  simulations  of  supersonic  and  hypersonic 
flows,  but  even  higher-order  methods  do  not  perform  correctly  in  turbulent  regions. 
Lee  et  al  [35]  found  that  the  use  of  a  6*^-order  ENO  scheme  was  too  diffusive  to 
properly  capture  turbulent  features  in  a  problem  of  shock  /  turbulence  interaction. 

Rather  than  trying  to  develop  a  single  scheme  that  performs  well  in  any  region 
of  a  supersonic  turbulent  flow,  a  hybrid  approach  has  been  developed  to  enable  such 
simulations.  The  base  solver  is  a  low-dissipation  central  scheme,  adapted  to  the 
simulation  of  turbulent,  reacting  flows,  which  has  extensively  and  successfully  been 
used  for  practical  studies  of  high  Reynolds  number  reacting  flows.  This  scheme  is 
combined  with  a  flux-difference  splitting  approach  that  permits  the  capture  of  flow 
discontinuities  locally. 
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3.2.1  Description  of  the  flux-difference  splitting  scheme 

The  upwind  scheme  used  in  the  present  formulation  is  a  high-order  flux  difference 
splitting  method.  The  reconstruction  of  the  interface  is  performed  using  a  MUSCL 
type  interpolation  [54]  with  a  limiting  function  to  ensure  monotonicity.  The  mono- 
tonized  central  limiter  is  adopted.  The  reconstruction  of  a  variable  U  consists  in 
extrapolating  the  variable  to  the  interface  with  a  directional  bias,  so  that  two  states 
are  obtained,  left  {U^)  and  right  {U^),  following  the  approach  given  by: 

^i+l/2  =  Ui-\ - 

Ui+l/2  ~  Ui+1  -  (j)mc{f'i^^/2)^i+3/2{U) 

where  Ai+i/2{U)  is  given  by: 


=  f/.+i  -  U, 


The  limiting  function,  used  here  is  the  monotonized  central  limiter,  which  satishes 
the  symmetry  condition  as  well  as  the  TVD  property  [27].  It  is  given  as: 


4>mc{r) 


max  [O,  min{2r,  2, 


1  -|-  r 
2 


)] 


where  is  a  coefficient  determined  by  a  flattening  operation  [12]  used  to  prevent  steep 
and  strong  gradients  from  creating  unphysical  oscillations  in  the  post-shock  region. 
This  flattening  is  done  by  redncing  the  order  of  reconstrnction  on  all  variables  when 
steep  strong  shocks  are  detected.  A  cell  is  identihed  as  being  part  of  a  shock  wave  if 
the  following  two  conditions  are  satished: 


\Pi+l  ~  Pi-l\  ^  1 

mm(Pi+i,Pi_i)  3 

(3.46) 

Ui+i  -  Ui_i  <  0 

(3.47) 

The  shock  thickness  is  then  measured  by  relating  the  pressure  gradient  across  two 
zones  to  the  gradient  across  fonr  zones. 


^  . 
^P,l 


Pi+l  —  Pi-1 
Pi+2  —  Pi-2 


(3.48) 


=  max 

The  limiting  factor  is  then  defined  as: 

max{^,  5+i) ,  ifPi+i  -  Pi-1  <  0 
max(^^i,^i-i),  otherwise 


0,  min (^l,10(^p,i- 0.75)) 


(3.49) 


(3.50) 
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The  reconstruction,  equations  3.46  and  3.46  is  then  fully  dehned,  and  can  be 
applied  to  various  sets  of  variables  ([4]).  In  the  present  work,  the  operation  is  per¬ 
formed  on  the  primitive  variables  (p,  Ui,  P,  pk)-  Additional  variables  if  needed  are  then 
computed,  and  the  monotonicity  of  these  derived  variables  is  explicitly  checked  and 
enforced.  The  monotonicity  of  the  interface  temperatures,  for  instance,  are  enforced 
by  reducing  the  order  of  extrapolation  of  pressure  and/or  density  as  needed.  The 
reconstruction  of  the  species  specihc  density  is  also  performed  with  special  attention. 
The  most  limiting  factor  of  all  specihc  densities  is  applied  to  all  others,  in  order  to 
ensure  mass  conservation  at  the  cell  interfaces. 

Fluxes  at  the  cells  interfaces  are  then  obtained  through  the  resolution  of  the  Rie- 
mann  problem  formed  by  the  reconstructed  state.  The  problem  is  solved  using  an 
approximate  Riemann  solver.  In  this  work,  a  hybrid  approximate  Riemann  solver 
has  been  developed,  that  permits  a  proper  capture  of  the  discontinuities  in  super¬ 
sonic  hows,  while  reducing  the  instabilities  of  other,  classical,  Riemann  solvers.  Its 
formulation  is  based  on  the  combination  of  two  solvers  of  the  HLL-  family  ([27]). 

The  HLLC  [53]  approximate  Riemann  solver  is  based  on  an  extension  of  the  2  — 
waves  HLL-Riemann  solvers  for  hydrodynamics  problems,  rather  than  on  the  formal 
3  —  waves  HLL—  solvers.  Its  formulation  results  in  a  computationally  inexpensive, 
robust  and  yet  accurate  solver  for  the  resolution  of  hydrodynamics  problems,  and 
compares  well  to  other  commonly  used  Riemann  solvers  (such  as  the  Roe,  Osher,  Two- 
Shock,  Two- Rarefaction  Riemann  solvers,  see  for  e.g.  [29],  [21])  It  is  often  preferred 
to  the  Roe-solver  for  its  characteristics  of  positivity  and  entropy  preservation.  The 
known  instability  of  this  solver  to  odd-even  decoupling  and  carbuncle  phenomenon  is 
cured  by  employing  the  HLLE  Riemann  solver  ([18]),  a  truly  2  —  waves  HLL-solver 
for  hux  computation  in  the  shock  transverse  directions  (detected  from  equations 
3.46  and  3.47).  This  hybrid  solver  has  a  reduced  sensitivity  to  instabilities,  and 
is  found  less  dissipative  than  hybrid  solvers,  which  use  HLLE  for  all  directions  in 
the  shock  ([45]).  In  particular,  the  new  solver  recovers  the  accuracy  of  the  HLLC 
Riemann  solvers  for  the  resolution  of  one-dimensional  problems.  Extension  of  the 
computational  method  to  multi-dimensional  problems  is  performed  using  an  unsplit 
methodology.  The  hybrid  Riemann  solver  does  not  show  any  instability  to  the  odd- 
even  decoupling,  and  accurately  resolves  the  carbuncle  test  case  [10].  More  details  on 
the  HLLE  and  HLLC  formulations  are  given  in  Appendix  A 

This  method  is  integrated  in  the  hybrid  formulation  by  sensing  strong  discontinu¬ 
ity  regions  using  a  smoothness  sensor.  The  switch  used  in  this  study  is  similar  to  that 
presented  in  [22]  and  applied  to  practical  turbulent  problems  in  [24].  The  smooth- 
flow  central  scheme  in  LESLIE3D  performs  well  for  the  resolution  of  turbulent  flows, 
but  cannot  correctly  simulate  strong  discontinuities.  A  smoothness  parameter  Si  is 
formulated  based  on  variables  curvatures  ([39]): 

o  _  \Qi+l  ~  ‘^Qi  T  Qi-l\  /„  N 

where  Q  can  be  any  variable  of  interest.  Both  pressure  and  density  can  be  used 
to  sense  the  smoothness  of  the  flow.  To  prevent  switching  on  numerical  noise,  a 
threshold  value  for  the  denominator  is  dehned.  The  smoothness  Si  is  set  to  0  if 
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either  the  numerator  or  the  denominator  is  less  than  O.OGQj.  The  cell  is  identihed 
as  non  —  smooth  if  the  smoothness  parameter  exceeds  a  threshold  value  of  0.5.  In 
multi-dimensional  problems,  the  transverse  fluxes  to  that  cell  are  also  tagged  as  non¬ 
smooth,  and  all  fluxes  are  evaluated  with  the  shock-capturing  scheme.  In  all  other 
cases,  the  smooth  spatial  discretization  is  used. 

3.2.2  Validation  of  the  flux-difference  splitting  scheme 

The  discontinuity-capturing  method  has  been  extensively  tested  and  validated  for  a 
wide  variety  of  flow  conhgurations.  Some  representative  tests  are  discussed  here. 

The  shock- capturing  capability  is  illustrated  by  iVoh’s  problem.  This  test  case 
simulates  two  infinite  strength  shocks  moving  outwards  from  the  center  of  a  one¬ 
dimensional  domain.  The  initial  conditions  are  given  by  a  sharp  separation  at  a;  =  0.5, 
the  center  of  the  domain,  where  the  left  state  is  given  by  (p,  m,  P)  =  (1, 1, 10“®),  while 
the  right  state  is  {p,u,P)  =  (1,— 1,10“®).  The  simulation  is  run  for  a  total  time 
T  =  1.  The  problem  is  solved  over  a  grid  with  100  grid  cells.  Figure  3.1  shows  the 
density  prohle  at  the  end  of  the  simulation  using  the  current  solver,  and  compared 
to  a  reference  solution.  In  the  course  of  the  simulation,  the  shocks  are  captured  at 
the  right  propagation  speed.  Furthermore,  the  post-shock  prohles  show  very  small 
oscillations,  as  a  result  of  the  flattening  operation,  leading  to  a  proper  resolution  and 
proper  shock  propagation  speed  with  limited  dissipation.  Only  a  slight  dip  is  observed 
in  the  density  prohle,  as  a  result  of  the  transience  in  the  formation  of  the  shocks. 


3 

a. 
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Figure  3.1:  Density  prohle  for  two  outwards-going  inhnite  strength  shocks. 


- 

- 

-  Reference  solution 
HLLC/E 

- 

_ ^ ^ _ , _ 

- 

The  extension  to  multidimensional  how  is  demonstrated  through  the  simulation 
of  a  point  explosion  problem  of  Sedov  on  a  two-dimensional  grid.  256  x  256  grid  cells 
are  used  to  discretize  a  Im  x  Im  domain  with  a  uniform  Cartesian  grid.  A  cylindrical 
region  of  overpressure  is  initialized  at  the  origin  (bottom  left  corner)  of  the  domain, 
over  a  radius  of  3.33cm.  The  pressure  ratio  across  the  initial  discontinuity  is  3.5  10®. 
The  strong  blast  wave  generated  by  this  energy  deposition  soon  reaches  a  state  of 
self-similar  evolution.  A  high-resolution  one-dimensional  simulation  of  this  problem 
in  cylindrical  coordinates  is  also  performed  for  reference.  Figure  3.2  represents  the 
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one-dimensional  projection  of  the  (normalized)  pressure  of  each  grid  cell  versus  the 
radius  of  the  cell  from  the  origin.  The  small  scattering  observed  in  this  plot  shows 
that  the  blast  wave  front  is  captured  at  its  right  speed,  even  when  its  propagation  is 
not  aligned  with  the  cells.  Furthermore,  the  capture  of  the  shock  front  is  relatively 
crisp,  extending  over  two  to  three  cells. 


Figure  3.2:  Sedov’s  point  explosion  pressure  prohles. 


The  HLLC/E  solver  described  earlier  has  been  designed  to  reduce  the  sensitive¬ 
ness  of  the  HLLC  to  a  weakness  shared  by  most  classical  Riemann  solvers.  The 
odd-even  decoupling  and  carbuncle  effect  are  solver’s  instabilities,  that  can  lead  to 
strong  amplitude  (unphysical)  oscillations  in  the  post-shock  regions,  which  strongly 
perturb  the  flow,  and  can  wrongly  interact  with  turbulent  structures.  The  Odd-Even 
decoupling  test  described  by  Liska  and  Wendroff  [37]  is  repeated,  using  the  HLLC/E 
solver.  Strictly  no  instability  is  detected  in  that  case,  whereas  similar  simulations 
using  the  HLLC,  the  Two-Shock  Riemann  Solver  and  the  Roe  Riemann  solver  with 
Harten-Hyman  entropy  corrections  showed  very  strong  sensitiveness  to  the  instability. 

The  classical  test  that  highlights  the  carbuncle  phenomenon  has  been  conducted: 
the  simulation  of  a  M  =  10  flow  over  a  circular  blunt  body  is  performed.  A  grid  with 
80x160  cells  is  used  to  resolve  this  problem.  Figure  3.3  shows  the  temperature  iso-lines 
from  the  computations  using  the  HLLE,  the  HLLC,  the  Roe,  the  Two-Shock  Riemann 
Solver  (TSRS)  [52],  and  the  HLLC/E  solvers.  The  result  obtained  with  the  HLLE 
solver  can  be  used  as  reference,  while  the  other  computations  show  an  instability  along 
the  main  shock  front,  leading  to  the  formation  of  spurious  oscillations  in  the  post¬ 
shock  region.  The  solver  of  Roe  is  particularly  sensitive  to  this  instability,  whereas 
the  HLLC/E  does  show  a  small  kink  within  the  shock  front,  but  its  instability  appears 
as  very  small  when  compared  to  the  other  solvers.  This  study  validates  the  HLLC/E 
hybrid  solver  and  its  reduced  instability. 

3.2.3  Validation  of  the  LES  methodology 

The  algorithms  developed  in  the  current  effort  and  described  earlier  have  been  used  to 
study  high-speed  compressible  turbulence  using  both  DNS  and  LES,  as  well  as  some 


18 


(a)  HLLE 


(b)  HLLC 


(c)  HLLC/E 


(d)  Roe 


(e)  TSRS 


Figure  3.3:  Carbuncle  test  case  using  contact-preserving  and  mixed  Riemann  Solvers 


canonical  test  problems.  The  hybrid  scheme  designed  to  capture  strong  discontinu¬ 
ities  and  to  resolve  the  smooth  variations  is  used  to  study  the  Shu  —  Osher  [51]  test 


problem.  This  test-case  consists  in  a  one-dimensional  shock  front  propagating  into  a 
sinusoidal  density  distribution.  The  interaction  forms  a  region  of  rapid  oscillations  in 
the  post-shock  region  which  steepen  into  shocks,  forming  an  N-wave  pattern.  A  cor¬ 
rect  resolution  of  this  problem  requires  an  accurate  computation  of  shock  propagation 
at  the  right  speed,  and  a  proper  capture  of  the  short-wavelength  smooth  variations 
that  form  in  the  post-shock  region. 

The  initial  conditions  for  this  problem  are  given  by  the  following  density  prohle: 


X  <  2 

X  >  2 


(3.52) 


The  traveling  shock  is  initialized  at  x  =  2,  and  the  pre-shock  pressure  and  velocity 
are  set  1  and  0,  respectively.  The  post-shock  values  for  pressure  and  velocity  are 
constant  and  equal  to  10.333333  and  2.629369,  respectively.  A  calorically  perfect  gas 
is  simulated  with  an  adiabatic  index  7  of  1.4.  The  simulation  is  run  for  a  total  time 
of  1.872. 

The  purely  upwind  approach  cannot  capture  short  wavelength  oscillations  when 
200  grid  cells  are  employed  (Fig.  3.4),  and  in  this  region  the  dominant  wavenumber 
of  the  N-wave  pattern  appears  as  the  smallest  resolved  wavenumber.  Employing  the 
hybrid  method  permits  for  the  same  resolution  to  crisply  capture  the  shock  with  the 
upwind  scheme,  while  the  smooth  flow  solver  is  used  to  resolve  the  post-shock  medium. 
As  a  result,  despite  the  fact  that  the  full  amplitude  of  the  oscillations  is  not  totally 
recovered,  the  short  wavelengths  are  apparent  in  the  post-shock  region.  Increasing 
the  resolution  to  400  grid  cells  conhrms  the  low-resolution  observation.  The  upwind 
approach  captures  all  the  wavelengths  of  the  problem,  but  the  inherent  dissipation 
of  the  scheme  prevents  this  approach  from  capturing  the  oscillations  amplitude.  The 


19 


hybrid  approach  combines  the  advantages  of  both  schemes,  leading  to  a  proper  capture 
of  the  main  shock  and  a  very  good  resolution  of  the  post-shock  oscillations. 


Figure  3.4:  Evolution  of  the  density  for  the  Shu-Osher  problem.  — ,  Hybrid 
(N=200) - ,  pure  upwind  (N=200)  ■  •  •,  Hybrid  (N=400)  o,  pure  upwind  (N=400)  o 


The  study  of  multi- mode  Richtmyer- Meshkov  Instabilities  is  shown  next.  A  do¬ 
main  of  0.04cm  X  0.01cm  is  discretized  using  1024  x  256  grid  cells.  An  initial  interface 
between  hydrogen  at  an  initial  temperature  of  SOOA'  and  helium  at  dhOOA'  is  shaped  by 
superposing  a  primary  perturbation  (of  wavelength  0.01cm  and  amplitude  0.001cm) 
and  a  secondary  perturbation  (with  a  wavelength  of  0.002cm  and  Xjlth  of  the  primary 
amplitude).  A  shock  with  an  overpressure  ratio  of  300  propagates  in  the  region  con¬ 
taining  helium,  towards  the  shaped  interface.  The  interaction  leads  to  a  180°  phase 
inversion.  The  initially  smooth  interface  forms  large  and  small  scale  vortices  under 
the  effect  of  baroclinic  torque,  which  transition  to  turbulence  and  eventually  decay. 
The  shock  is  partially  transmitted  in  the  hydrogen  region,  and  is  regularly  distorted 
by  a  series  of  transverse  waves  that  form  local  triple  points.  The  slip  lines  so  formed 
form  a  diamond  —  shaped  pattern,  which  creates  vortical  structures  due  to  Kelvin- 
Helmholtz  instabilities.  This  problem  has  been  solved  using  the  hybrid  methodology 
with  both  the  Piecewise  Parabolic  Method  (PPM)  and  with  the  current  HLLC/E  up¬ 
wind  scheme.  The  PPM  implementation  is  identical  to  that  of  the  PROMETHEUS 
code  and  the  FLASH  code  developed  in  the  university  of  Chicago  [23],  and  its  in¬ 
tegration  in  the  hybrid  methodology  has  been  performed  and  described  earlier  [22]. 
Fig.  3.5  shows  the  temperature  helds  for  all  computations,  and  the  regions  where  up- 
winding  is  used  in  the  fluxes  computation.  The  large  structures  formed  through  the 
interaction  are  found  identical  with  both  methods.  The  small  scale  vortices  are  found 
to  be  slightly  different  for  the  two  cases,  but  comparable.  The  regions  of  upwinding 
are  signihcantly  more  important  in  the  case  when  PPM  is  used  for  discontinuity  cap¬ 
turing  purposes,  as  illustrated  in  Fig.  3.5.  The  directional  splitting  used  in  the  PPM 
method  requires  an  updated  of  all  the  neighboring  cells,  hence  broadening  the  stencil 
of  PPM  usage  in  the  hybrid  formulation.  The  unsplit  HLLC/E  method,  on  the  other 
hand,  used  the  same  integration  method  as  the  smooth  flows  solver,  and  can  then  be 
locally  activated.  It  leads  then  to  a  smaller  region  of  activation  around  the  density 
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interfaces,  hence,  a  reduced  use  of  the  dissipative  scheme. 


Hybrid  /  PPM 


Hybrid  /  HLLC/E 


Hybrid  /  PPM 


. 

Hybrid  /  HLLC/E 


(a)  (b) 

Figure  3.5:  Simulations  of  multi-mode  Richtmyer- Meshkov  instabilities  using  hybrid 
PPM/4*^  (top)  and  hybrid  HLLCE/4*^  (bottom):  Temperature  helds  (left)  and  switch 
helds  (right). 


The  Richtmyer-Meshkov  Instability  (RMI)  problem,  used  here  to  highlight  the 
switching  process  and  for  comparison  of  the  hybrid  methods,  is  a  fundamental  insta¬ 
bility  in  shock  containing  flows.  The  mixing  induced  through  the  interaction  and  the 
turbulence  generated  by  the  long-time  decay,  and/or  by  a  re-shock  make  this  case  a 
computational  challenge.  This  conhguration  will  be  re-examined  later  using  LES  as 
well. 

The  hybrid  scheme  can  reliably  capture  shocks  and  still  apply  a  limited  amount 
of  dissipation  on  the  turbulent  structures.  The  closure  model  is  evaluated  by  doing 
LES  and  comparing  with  DNS  using  the  current  code  and  also  with  reference  DNS 
results  [40],  [35].  The  DNS  studies  have  been  used  in  a  priori  studies  to  highlight 
the  driving  terms  of  the  governing  equation  within  this  conhguration.  The  LES 
grid  used  to  discretize  the  conhguration  shown  in  Fig.  3.6  is  coarser  by  a  factor  of  13 
compared  to  the  DNS  grids.  The  resolution  in  the  direction  across  the  shock  region 
is  the  same  spacing  as  in  the  DNS  in  order  to  capture  the  shock  but  the  resolution 
in  the  other  two  directions  were  reduced  drastically.  The  initializations  of  the  LES 
simulation  is  obtained  from  the  hltered  DNS  held. 

The  time-averaged  axial  Reynolds  stress  prohles  obtained  from  the  calculations  of 
the  M  =  1.29  and  M  =  3.0  cases  are  shown  in  Fig.  3.7  along  with  their  corresponding 
hltered  DNS  helds.  In  these  plots,  /cq  denotes  the  most  energetic  wavenumber  in  the 
initial  turbulent  held  used  at  the  inhow,  and  is  used  to  normalize  the  distance  from 
the  inhow.  The  results  of  the  LES  using  the  LDKM  approach  are  shown  along  with 
the  results  of  an  under-resolved  simulation  without  the  closure  model.  The  turbulent 
decay  that  occurs  before  the  interaction  is  well  captured  at  the  right  rate,  when  the 
explicit  subgrid  model  is  used.  Without  model,  the  turbulent  decay  is  oh  for  the  low 
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Mach  number  case.  No  real  trend  can  be  isolated  from  the  highest  Mach  number  case, 
as  the  timescale  for  convection  of  the  turbulent  structures  is  significantly  smaller  than 
the  timescales  for  production  and/or  dissipation.  The  post-shock  level  of  stresses  is 
correctly  predicted  by  the  numerical  method.  The  peak  values  reached  downstream 
the  shock,  from  the  acoustic  /  kinetic  energy  transfer,  are  slightly  under-estimated  in 
both  cases.  Even  though  the  minimum  spacing  matches  that  of  the  DNS  simulation, 
the  grid  around  the  shock  is  coarser,  and  leads  to  a  slight  under-prediction.  The 
dynamic  model  is,  however  seen  to  reproduce  the  gain  in  turbulence  in  the  post¬ 
shock  region  correctly.  The  turbulent  decay  that  follows  is  captured  at  the  right 
rate  when  the  LDKM  is  used.  The  absence  of  model,  on  the  other  hand,  fails  to 
capture  the  turbulent  behavior.  The  molecular  viscosity  in  this  under-resolved  does 
not  provide  enough  dissipation  to  reproduce  the  hltered  DNS  prohles. 


Figure  3.6:  Computational  set-up  for  the  simulation  of  shock  /  isotropic  turbulence 
interaction.  A  box  of  isotropic  turbulence  is  used  at  the  inflow  of  a  spatial,  shock-hxed 
domain. 


Figure  3.7:  Longitudinal  Re-stress  prohles  for  a  M  =  1.29  Shock  (left)  and  a 
M  =  3.0  shock  (right)  interacting  with  an  isotropic  turbulent  held. 
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The  results  presented  in  this  section  highlight  the  different  aspects  of  the  algo¬ 
rithms  developed  in  this  research.  The  upwind  scheme  permits  a  proper  capture  of 
the  shock  waves  over  2/3  cells,  and  with  the  proper  propagation  speed.  The  hy¬ 
brid  scheme  combines  this  upwind  method  to  a  smooth  central  scheme  adapted  to 
turbulent  flows  simulation.  This  hybrid  scheme  successfully  captures  the  problem 
of  Shu-Osher,  shock  /  turbulence  interaction.  Furthermore,  the  development  of  the 
LDKM  closure  approach  have  permitted  a  good  capture  of  the  turbulent  behaviors  for 
this  interaction  without  requiring  any  ad  hoc  model  adjustment.  In  fact,  the  entire 
algorithm  (both  numerical  and  subgrid)  developed  in  this  effort  requires  no  model 
adjustments  in  any  of  the  reported  flows  (see  next  section).  The  application  of  this 
method  to  different  flows  of  practical  interest  is  presented  next. 
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CHAPTER  IV 


RESULTS  AND  DISCUSSION 


4.1  Richtmyer- Meshkov  Instability 

A  Richtmyer-Meshkov  Instability  (RMI)  is  the  result  of  an  impulsive  acceleration  of 
an  interface  between  two  media  of  different  densities.  The  initial  density  discontinuity 
leads  to  the  formation  of  a  momentum  interface  and  of  a  mixing  region  which  grows 
linearly  in  the  early  time  of  the  experiment.  When  non-linearities  start  playing  a 
major  role,  the  large-scale  structures  start  breaking  down,  and  a  region  of  turbulence 
is  formed  from  the  initial  interface. 

Classically,  the  impulse  needed  in  experimental  studies  of  Richtmyer-Meshkov 
instabilities  is  provided  by  a  traveling  shock.  As  the  shocks  travel  down  the  shock 
tubes,  they  reflect  from  the  back  wall  of  the  experimental  facility,  and  re  —  shock  the 
mixing  region.  This  was  observed  to  accelerate  the  transition  to  turbulence  for  this 
instability.  One  of  the  experiments  presented  by  Vetter  and  Sturtevant  [55]  has  been 
studied  numerically  with  the  current  LES  methodology.  A  M  =  1.98  shock  interacts 
with  an  Air/SF^  interface,  leading  to  an  initial  RMI,  which  forms  large  coherent 
structures.  The  shock  reflection  from  the  back-wall  re-shocks  the  mixing  zone,  and 
transition  to  turbulence  is  enhanced.  Experimental  measurements  of  the  mixing  zone 
are  used  for  comparison.  The  domain  employed  for  the  current  simulation  is  a  0.59m 
long  with  a  0.27m^  cross  section  shock  tube,  and  is  schematically  shown  in  Fig.  4.1. 
Periodic  conditions  are  enforced  in  y  and  z  direction.  The  right  boundary  was  set  as 
a  no-flux  wall,  whereas  the  left  boundary  is  treated  as  an  outflow. 


Figure  4.1:  Geometry  of  the  Simulation 

In  the  numerical  simulation,  the  interface  between  air  and  SFq  is  initialized  as: 
Xi{y,z)  =  Oo  |sin(7r|//A)  sin(7r2:/A)|  +  aih{y,  z),  where  h(y,z)  is  a  random  function, 
which  perturbs  the  initial  interface  prohle  to  break  the  symmetry  and  enhance  per¬ 
turbation  and  transition.  Here,  oq  =  0.25cm,  ai  =  0.025cm  and  A  =  0.27/14cm  are 
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Case 

1 

2 

3 

4 

Physical  Size  (m) 

0.59x0.27^ 

0.59x0.27^ 

0.59x0.272 

0.59x0.1352 

Grid  Size 

281x129^ 

420x1932 

560x2572 

560x1292 

Grid  Resolution  (cm) 

0.21 

0.14 

0.105 

0.105 

Simulation  Time  (ms) 

8.0 

6.0 

6.0 

6.0 

CPU  hours 

780 

1680 

7360 

2040 

Table  4.1:  Grid  resolution  and  domain  size  for  various  test  cases 


Figure  4.2:  Mixing  length  before  reshock 


Figure  4.3:  Overall  mixing  length  growth 


assumed  [28].  Resolution  and  domain  size  studies  have  been  performed.  Four  test 
cases  (Tabled.  1)  are  discussed  here. 

Vetter  and  Sturtevant  [55]  estimated  the  instantaneous  mixing  lengths  from  visual 
estimates  of  Schlieren  images.  The  numerical  mixing  length  is  defined  as  h{t)  = 
4  <  'ip  >)  <  'ip  >  dx,  where  <>  represents  the  planar  average.  The  growth 

rate  before  reshock  can  be  often  represented  as  hi(t)  =  ar^  , where  h  represents  the 
mixing  length,  a  is  a  constant,  and  r  is  linear  function  of  the  time  [59].  Zhou  [60] 
found  that  db  varies  between  0.2  and  2/3.  The  present  simulation  falls  in  the  range  of 
predicted  values,  with  6  =  0.59  (Fig.  4.2).  Figure  4.3  clearly  shows  that  the  mixing 
width  is  over-predicted  by  the  coarse  grid  but  Cases  3  &  4  show  that  the  mixing 
length  do  not  depend  on  the  largest  scale  of  physical  domain,  and  that  a  reduced 
physical  length  can  be  used.  There  is  good  agreement  with  the  experimental  data 
|66]. 

A  time  sequence  of  iso-vorticity  contours  is  shown  in  Fig.  4.4.  Figs.  4.4(a)  and 
4.4(b)  show  the  structural  evolution  from  the  first  shock  /  interface  interaction.  The 
formation  of  large  finger— like  shapes  is  initiated  from  the  instability.  These  struc¬ 
ture  grow  in  the  linear  regime,  and  very  high  levels  of  vorticity  are  obtained.  These 
structures  would  eventually  transition  to  turbulence  and  decay.  The  re-shock,  how¬ 
ever,  enhances  this  instability,  and  Figs.  4.4(c)  and  4.4(d)  show  the  vortical  structures 
after  the  second  passage  of  the  shock.  It  is  observed  that  the  large,  organized  struc¬ 
tures  have  broken  down  as  a  result  of  the  interaction,  forming  a  wide  mixing  region 
with  numerous  vortical  structures  of  different  scales.  2D  spectra  of  the  axial  velocity 
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(c)  4.0ms,  8000(l/s) 


(d)  5.0ms,  8000(l/s) 


Figure  4.4:  Iso-surface  of  total  vorticity  magnitude,  colored  by  mass  ratio  of  air 


26 


(a)  Before  reshock  (t=1.5ms)  (b)  After  reshock  (t=5.0ms) 

Figure  4.5:  2D  Spectra  of  the  x-component  of  velocity 


within  the  mixing  zone  (dehned  as  the  region  where  V’  >  0.1)  are  presented  in  Fig. 
4.5.  The  dominant  wavenumber  imposed  at  the  initialization  and  its  harmonics  are 
persistent  in  the  spectra  shown  in  Figs.  4.5(a).  This  is  in  accordance  with  the  large- 
scale  structures  observed  from  the  iso-vorticity  contours.  Again,  the  late  development 
would  dissipate  the  dominant  waves  and  allow  the  formation  of  turbulent  structures. 
However,  re-shocking  is  achieved  before  this  transition  stage  is  reached.  Figure  4.5(b) 
shows  the  spectrum  obtained  after  re-shock.  The  second  interaction  enhances  the 
mixing,  and  a  wider  spectrum  is  obtained  without  dominant  frequency.  Zhou  [59]  an¬ 
alytically  derived  that  the  energy  spectrum  of  the  RMI  shows  a  scaling,  rather 
than  the  classical  due  to  inhomogeneity  and  anisotropy  induced  by  the  shock 

motion.  In  the  spectra,  an  inertial  range  is  obtained,  although  the  dominant  slope 
does  not  match  exactly  the  theoretical  values  or  but  is  closer  to 

This  latter  scaling  was  also  obtained  in  another  LES  study  [11]  but  is  currently  still 
being  investigated. 

The  RMI  is  a  challenging  problem  in  numerical  studies,  as  it  evolves  from  a  sharp 
density  interface.  As  such,  the  scheme  is  required  to  capture  strong  discontinuities 
with  limited  dispersion,  and  still  show  small  dissipation  in  order  to  capture  the  most 
unstable  modes  of  the  conhguration.  The  late-stage  developments  of  instability  lead 
to  the  formation  of  high  strain  and  high  vorticity  regions,  which  eventually  form 
turbulence.  As  noted  before,  re-shocking  the  interface  only  accelerates  this  process. 
This  second  stage  requires  a  high  resolution  scheme  with  low  dissipation,  and  a  correct 
subgrid  modeling  that  will  mimic  turbulence  generation  and  the  subsequent  turbulent 
decay. 

4.2  Regular  and  Irregular  Cellular  Detonations 

In  a  detonation  wave,  a  high-speed  shock  travels  in  a  low-temperature  reactive  mix¬ 
ture,  compressing  and  heating  the  gases.  In  the  post-shock  region,  the  conditions  to 
achieve  auto-ignition  are  reached.  The  heat  released  by  the  reaction  and  the  gas  ex¬ 
pansion  feed  back  into  the  main  wave.  A  good  capture  of  the  phenomenon  requires  a 
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proper  resolution  of  the  running  shock,  an  accurate  treatment  of  the  chemical  source 
terms,  and  a  low  dissipation  resolution  in  the  post-shock  region  in  order  to  capture 
the  turbulent  structures  emanating  from  the  triple  points  on  the  detonation  front. 

Simulation  of  regular  and  irregular  cellular  detonations  is  challenging  as  it  involves 
large  range  of  length  and  time  scales,  nonlinear  coupling  between  chemistry  and  fluid 
dynamics  and  instability  in  flow.  Different  detonations  have  been  simulated,  and  two 
cases  are  presented  hereafter,  one  of  regular  and  one  of  irregular  nature  (see  Table  4.2). 
The  initial  prohles  of  pressure,  temperature  and  species  mass  fractions  are  obtained 
using  a  one-dimensional  steady-state  (ZND)  detonation  solver.  These  prohles  are 
projected  on  a  two-dimensional  grid,  and  the  perturbation  that  transitions  to  the 
cellular  structure  formation  is  initiated  by  placing  a  pocket  of  unreacted  mixture 
behind  the  front  ([43]). 


Case 

Type 

Resolution 

Cell  resolution 
{Ax  =  Ay)  fj,m. 

CPU  hours  used 

I 

Regular 

4800X1200 

7.25 

1665 

II 

Irregular 

4000X2000 

18.75 

2850 

Table  4.2:  Resolution  for  regular  and  irregular  detonation  simulations. 

First,  a  detonation  with  regular  cellular  structures  is  presented:  Case  I  considers 
the  regular  detonation  obtained  far  a  stoichiometric  mixture  where  7  is  1.2,  the  initial 
temperature  (Tq)  293K,  the  initial  pressure  0.2atm.,  the  activation  energy  (Ea/RTo) 
is  50  and  the  heat  released  by  the  reaction  is  given  by  Q/RTq  =  50.  Here,  R  is 
universal  gas  constant.  The  Chapman-Jouguet  (CJ)  detonation  speed  {Dcj)  is  the 
most  likely  detonation  speed,  and  is,  for  this  case,  1.8km/s.  The  overdrive  parameter 
(/  =  {D/DcjY)  is  1.6,  where  D  is  actual  detonation  speed.  Recording  the  spatial 
evolution  of  the  triple  points  tracks  leads  to  a  numerical  sootfoil  track.  Figure  4.6 
shows  the  regular  cells  formed  for  Case  I.  This  track  is  compared  to  the  numerical 
results  of  Deiterding  [14]  for  the  same  conditions. 


(a)  (b) 


Figure  4.6:  Numerical  soot  foil  recording  maximum  pressure  showing  regular  cellular 
structures  from  the  current  simulation  (left)  and  from  Deiterding  [14]  (right). 


Case  II  is  for  an  irregular  detonation.  The  physical  parameters  for  this  case 
are:  7  =  1.228,  Ea/RTo  =  38.2  and  Q/RTq  =  34,  and  the  overdrive  factor  is  /  = 


1.4.  Figure  4.7  shows  the  sootfoil  for  this  conhguration  along  with  an  instantaneous 
capture  of  the  detonation  front,  represented  through  a  snapshot  of  the  held  of  density 
gradient.  The  sootfoil  shows  the  irregularity  of  the  cells  formed  by  the  triple  points. 
Some  triple  points  fade  out  as  a  consequence  of  the  weak  transverse  waves  induced  in 
this  case.  The  instantaneons  held  also  shows  some  of  the  vortical  strnctnres  indnced 
by  the  triple  point  from  baroclinic  ehects  and  Kelvin-Helmholtz  instabilities.  The 
complex  nature  of  the  irregular  detonation  comes  from  the  interaction  between  the 
transverse  waves,  the  vorticity  and  the  heat  released  from  the  reaction. 


Figure  4.7:  Nnmerical  soot  foil  (left)  and  instantaneons  held  of  density  gradient 
(right)  for  an  irregnlar  detonation  (case  II) 


Fnll  3D  LES  of  detonation  has  also  been  stndied.  A  case  of  detonation  is  simulated 
in  a  mixture  of  2H2/102/^Ar  initially  at  a  pressnre  of  2QkPa  and  a  temperatnre  of 
30077,  similar  to  the  experiments  of  [6].  A  7-steps  H2/O2  reaction  mechanism  is 
employed  to  treat  the  combustion  [20],  and  500  x  180  x  180  grid  cells  are  used  to 
discretize  a  physical  domain  of  dimensions  8cm  x  1.25cm  x  1.25cm.  Fignre  4.8  shows 
a  snapshot  of  the  detonation  strnctnre.  The  transverse  waves  in  two  dimensions 
are  clearly  identihed  from  the  pressnre  gradient  held.  Their  regnlar  strnctnre  in  three 
dimensions  is  illustrated  by  the  regular  pattern  of  instabilities  that  corrugate  the  main 
front  of  iso-pressnre.  Also,  the  nnsteady  triple  points  create  vortical  strnctnres  in  the 
post-detonation  region  that  lead,  throngh  Kelvin-Helmholtz  instabilities,  to  tnrbnlent 
strnctnres  fnrther  downstream.  The  cell  width  obtained  from  this  simnlations  is 
predicted  in  close  agreement  to  the  experimental  observations. 

The  evolntion  of  the  detonation  front  in  three-dimensional  simnlations  reqnires  the 
nse  of  an  npwind,  shock-capturing  scheme.  Clearly,  the  downstream  evolntion  which 
depends  on  the  regnlarity  of  the  detonation  is  showing  nnstable  vortical  strnctnres, 
and  as  the  detonation  becomes  more  nnstable,  the  level  of  post-detonation  tnrbulence 
can  signihcantly  increase.  The  hybris  method  appears  well  adapted  to  captnre  of  snch 
complex  detonative  processes. 

4.3  Shock- Shear  Interactions  in  a  Scramjet 

Facing  the  problem  of  rednced  mixing  in  high-speed  shear  layers,  and  motivated  by  the 
experimental  and  analytical  prediction  of  tnrbulence  enhancement  throngh  shocks. 
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Figure  4.8:  Three-dimensional  LES  results  of  the  pressure  iso-surface  of  TOOfcPa, 
colored  by  the  temperature  held.  The  temperature  held  is  shown  in  the  last  plane, 
and  the  pressure  gradient  magnitude  is  on  the  bottom  plane. 


many  researchers  have  studied  the  development  of  shocked  turbulent  shear  layers. 
The  physical  processes  occurring  during  a  shock  /  turbulent  shear  layer  interaction 
are,  however,  numerous  and  not  limited  to  shock  /  turbulence.  The  vorticity  of  the 
large  coherent  structures  is  amplihed  through  the  strong  compression.  Also,  density 
gradients  across  the  shear  layer  can  lead  the  formation  of  local  Richtmyer-Meshkov 
instabilities  if  misaligned  with  the  shock’s  pressure  gradients.  Finally,  the  mean 
vorticity  within  a  shear  layer  does  not  necessarily  increase  through  the  interaction. 

The  conhguration  studied  numerically  is  that  of  a  two-streams  binary  mixing 
layer  ([15]).  The  upper  stream  (hereafter,  denoted  with  a  subscript  u)  is  a  mixture 
of  N2/H2,  with  10%  hydrogen  in  mass,  with  a  static  temperature  of  20007^  and 
a  static  pressure  of  latm,  flowing  at  Mach  2.0  (which  corresponds  to  a  velocity  of 
2672m /s  with  the  thermally  perfect  gas  Equation  Of  State  (EOS)  employed  here). 
The  lower  stream  (hereafter,  denoted  with  a  subscript  /)  is  that  of  air  where  the 
static  temperature  and  pressure  are  set  to  match  that  of  the  upper  stream.  The 
Mach  number  is  set  to  2.0,  which  corresponds  to  a  velocity  of  1729m/s.  Again, 
a  thermally  perfect  gas  EOS  is  used  for  this  stream  species.  A  schematic  of  the 
conhguration  studied  here  is  shown  in  Fig.  4.9  The  convective  velocity  for  this  how 
is  about  Uc  =  2100m/s.  The  convective  Mach  number  is  Me  =  0.43,  which  makes 
it  moderately  compressible,  with  turbulent  structures  that  still  resemble  those  of  the 
incompressible  mixing  layer.  In  order  to  trigger  the  transition  of  this  initially  laminar 
prohle,  a  random  phase  velocity  perturbation  is  added  to  the  mean  prohle  with  a 
hxed  energy  spectrum. 

The  problem  was  discretized  using  250  x  80  x  40  grid  cells  to  resolve  a  domain 
of  extent  0.17m  x  0.10m  x  0.03m.  Fig.  4.10  shows  a  top  view  of  the  iso-surfaces 
of  the  Q— criterion  for  the  base  mixing  layer  and  the  shocked  counter-part.  This 
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Figure  4.9:  Sketch  of  the  shocked  mixing  layer  conhguration. 


variable,  dehned  as  the  second  invariant  of  the  velocity  gradient  tensor,  is  well-suited 
for  vortical  held  identihcation  ([16]).  These  snapshots  are  taken  at  the  same  physical 
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Figure  4.10:  Iso-surface  of  the  Q-criterion  {Q  =  for  the  un-shocked  (top) 

and  the  shocked  (bottom)  mixing  layer,  colored  by  the  local  Mach  number  (how  is 
from  left  to  right) 


time,  after  10  how-through-times.  A  fast  transition  to  turbulence  is  initiated  by  the 
forcing  at  the  inhow  of  the  simulation,  leading  to  a  fully  developed  turbulent  mixing 
layer.  Strongly  three-dimensional  coherent  structures  are  observed  in  the  mixing 
layer.  These  coherent  structures,  initiated  by  natural  instabilities  of  the  how,  trigger 
and  maintain  the  turbulence  within  the  mixing  layer  thickness. 

The  shock  impact  is  found  to  amplify  all  the  turbulent  levels  within  the  mixing 
zone.  In  particular,  the  turbulent  huctuations  in  the  cross-wise  direction  and  the 
principal  Reynolds  stress  Rxy  are  signihcantly  enhanced  through  the  interaction.  As 
a  consequence,  the  vortical  structures  in  the  shocked  mixing  layer  are  compressed,  and 
show  more  two-dimensionality  than  in  the  un-shocked  case,  as  is  visible  in  Fig.  4.10. 
Later  downstream,  these  structures  recover  a  strong  three-dimensionality.  Further¬ 
more,  these  increases  in  turbulent  levels  lead  to  an  important  increase  in  the  mixing 
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efficiency  in  the  region  that  follows  the  shocks.  The  thickness  based  on  the  90%-iir2 
mass  fraction  is  shown  in  Fig.  4.11(a)  for  the  base  mixing  layer  and  for  the  shocked 
shear  layer.  A  reduction  of  the  thickness  is  observed  first  as  the  shock  interacts  with 
the  mixing  region.  This  is  due  to  the  spatial  compression  of  the  mixing  region  by  the 
two  shocks,  and  should  not  be  interpreted  as  a  reduction  in  mixing  efficiency.  On  the 
contrary,  the  shocked  layer  evolves  with  an  increased  growth  rate  in  the  region  that 
follows  the  interaction.  The  growth  rate  is  found  to  be  double  the  growth  rate  of  the 
un-perturbed  layer.  This  observation  is  conhrmed  by  the  prohle  of  mass  entrained  by 
the  mixing  layer,  showed  in  4.11(b).  The  mixing  achieved  through  the  shock  impact 
is  signihcantly  greater  than  in  the  case  of  the  base  mixing  layer. 


layer 

Figure  4.11:  Mixing  Layers  growth  rate  based  on  a  90%  H2  mass  fraction,  and  mass 
entrainment  within  the  mixing  layer  thickness 


Within  6cm  from  the  location  of  the  interaction,  the  mixing  layer  growth  rate 
steadies  out  at  a  value  that  is  about  6%  higher  than  its  unshocked  counter-part,  and 
the  turbulent  fluctuations  relax  to  their  un-disturbed  levels.  The  reduction  in  the 
mean  levels  of  vorticity  in  this  conhguration  leads  to  a  lower  amount  of  turbulent 
production  across  the  layer. 

The  effects  of  the  shocks  on  the  turbulent  evolution  and  on  the  mixing  efficiency 
are  found  to  be  strong,  but  quite  localized  in  space.  These  results  are  in  agreement 
with  some  experimental  observations  [50]  where  an  increase  in  the  growth  rate  was  ob¬ 
served  following  the  interaction  with  the  shocks,  with  a  relaxation  to  the  un-disturbed 
levels  observed  within  some  short  distance  from  the  point  of  impact.  It  should  also 
be  noted,  as  pointed  out  by  Drummond  et  at  [15],  that  the  shocks  interaction  also 
impact  the  static  temperature  within  the  layer,  hence  enhancing  the  reactions  rates 
in  the  mixing  zone.  Considering  that  shock  waves  cannot  be  avoided  in  practical 
Scramjet  engines,  this  method  can  be  considered  an  interesting  mixing-augmenting 
technique. 
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CHAPTER  V 


CONCLUSIONS  AND  FUTURE  PLANS 


A  LES  methodology  has  been  presented  that  permits  the  resolution  of  compressible 
turbulent  flows  in  supersonic  environments.  Two  key  developments  were  required 
to  achieve  this  goal:  a  new  hybrid  numerical  algorithm  that  captures  strong  shocks 
and  shear  turbulence,  and  a  localized  dynamic  subgrid  closure  for  compressible  flows. 
The  hybrid  numerical  algorithm  combined  a  shock  capturing  approach  with  a  shock- 
free  algorithm  using  an  automatic  shock  detection  technique  and  has  been  validated 
for  various  canonical  flows  containing  strong  shocks  as  wells  in  flows  with  shock- 
turbulence  interactions.  It  is  shown  that  this  algorithm  has  the  requisite  accuracy  to 
do  DNS  of  such  flows.  Subsequently,  an  advanced  closure  model  for  the  subgrid  kinetic 
energy  equation  has  been  developed  and  validated.  The  LDKM  approach  allows 
localized  evaluation  of  all  model  coefficients  during  the  simulation  and  therefore, 
there  are  no  adjustable  parameters  in  the  proposed  simulation  strategy. 

This  new  LES  approach  has  been  used  for  the  study  of  various  applications.  A 
case  of  re-shocked  Richtmyer-Meshkov  Instability  is  discussed  that  highlights  the 
importance  of  capturing  the  pressure  and  the  density  discontinuities  until  a  high 
level  of  vorticity  is  achieved  and  transition  to  turbulence  occurs.  Also,  fundamental 
studies  of  regular  and  irregular  detonations  are  conducted  showing  how  the  post¬ 
detonation  region  is  subject  to  the  formation  of  vortical  structures,  which  in  3D  lead 
to  the  formation  of  a  turbulent  medium.  Finally,  the  evolution  of  a  turbulent  shear 
layer  after  its  interaction  with  two  symmetric  shocks  is  studied.  An  increased  level 
of  turbulence  was  captured  behind  the  interaction  point,  leading  to  an  enhanced 
growth  rate.  This  phenomenon,  associated  to  the  increase  in  static  temperature 
make  this  mixing  augmentation  technique  apparently  well  suited  for  flame  intiation 
and  stabilization  in  supersonic  flows.  The  amplihed  turbulence  is,  however,  found  to 
decay  within  a  short  distance  from  the  interaction. 

The  LES  approach  developed  here  appears  well  adapted  to  the  simulations  of 
high-speed  turbulent  flows.  Its  applicability  to  practical  conhgurations  has  been 
demonstrated.  The  resolution  requirements  can,  however,  increase  signifcantly  as 
the  domains  become  more  and  more  complex.  A  next  step  in  the  developments  of 
this  LES  strategy  will  consist  in  implementing  an  Adaptive  Mesh  Rehnement  (AMR) 
technique.  This  method,  which  permits  the  grid  to  automatically  adjust  its  levels  of 
rehnement  is  ideally  suited  within  the  hybrid  framework  developed  in  this  effort. 
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APPENDIX  A 


HLLE  AND  HLLC  RIEMANN  SOLVERS 


A  brief  description  of  the  HLL  Riemann  solver  family  is  given  hereafter.  The  HLL 
approximate  Riemann  solvers  have  hrst  been  developed  by  Harten,  Lax  and  van  Leer 
([27]),  by  expressing  a  hyperbolic  system  of  conservation  laws  in  integrand  form.  From 
an  initial  interface  separating  two  constant  states,  it  is  assumed  that  N  waves  can  be 
formed  from  the  characteristics  evolution  of  the  system,  separating  iV  +  1  constant- 
properties  regions.  The  knowledge  of  the  jump  relation  through  the  waves  and  the 
wave-speeds  can  be  used  to  obtain  a  closed  form  expression  for  the  intermediate  states, 
and  the  associated  fluxes.  Harten  [27]  carried  out  a  full  derivation  for  a  2  —  waves 
problem,  leaving  the  wavespeeds  as  sole  unknowns  to  the  solver,  and  have  given  a 
mathematical  description  for  the  3-waves  problem. 

The  2-waves  formulation  for  the  resolution  of  the  Euler  equations  (with  wavespeeds 
expressions  given  by  [18,  19],  thus  called  HLLE)  has  been  proven  robust  and  adequate 
for  shocks  and  rarefactions,  but  appears  to  be  very  dissipative  for  contact  discontinu¬ 
ities.  In  another  study,  a  correction  to  the  derivation  of  this  solver  was  proposed  ([53]), 
in  order  to  add  the  missing  contact  wave  (thus,  called  HLLC),  whose  wavespeed  was 
estimated  by  an  approximation  of  the  particle  velocity  in  the  intermediate  region. 
The  formulation  is  closed  by  expressing  the  jump  conditions  across  all  wave  obtained 
from  the  exact  Riemann  solver  for  the  Euler  equations.  This  formulation  was  further 
studied  by  Batten  et  al  [3],  where  the  averaged  intermediate  state  was  related  to  the 
HLLE  evaluation,  thus  leading  to  an  easy,  but  yet  robust,  3- waves  Riemann  solver.  It 
should  be  noted  that  this  3-waves  solver  does  not  follow  the  original  work  of  Harten 
et  al,  but  rather  modihes  the  2-wave  system  by  introducing  the  intermediate  wave’s 
speed  from  the  2-wave  solver  as  a  correction.  The  solver  does  not  reduce  to  a  single¬ 
wave  problem  in  the  physical  limit  of  an  isolated  discontinuity.  It  is,  however,  found 
robust  and  yet  accurate  for  the  resolution  of  hydrodynamics  problem,  and  remains 
very  inexpensive  compared  to  other  approximate  Riemann  solvers. 

The  Riemann  solver  developed  for  this  study  uses  a  combination  of  both  solvers 
in  order  to  reduce  the  instabilities  associated  to  contact-resolving  solvers.  The  simple 
2-waves  HLL  will  be  described  hereafter.  The  wave-speeds  estimates  for  the  HLLE 
approach,  and  the  modihcations  that  lead  to  the  HLLC  solver  are  also  described. 

A.l  The  2- waves  HLL  Riemann  Solver 

This  approximate  Riemann  solver  starts  with  the  assumption  that  an  initial  dis¬ 
continuity  gives  rise  to  2  waves,  a  left-moving  wave  (of  Eulerian  speed  Si),  and  a 
right-moving  counterpart  (of  Eulerian  speed  Sr).  A  typical  (x,t)  diagram  for  a  sub¬ 
sonic  case  is  given  in  Fig.  A.l.  The  integral  form  of  the  Euler  equations  (see,  eg. 
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Figure  A.l:  Two  characteristic  waves,  and  their  characteristic  wavespeeds. 


[52]),  reads: 


F{U)dt]  =  0 


(A.l) 


where,  noting  q  =  V  ■  n  =  un^  +  vUy  +  wn^  the  interface  normal  velocity  amplitude, 
U  and  F{U)  are  given  by: 


p 

pq 

pu 

puq  +  Prix 

pv 

fyvq  +  Puy 

pw 

,  F  = 

pwq  +  Friz 

pE 

{pE  +  P)q 

pk^gs 

pk^gsq 

.  P^k  . 

pYkq 

Let  T  be  the  local  time  step,  T  =  >  0.  Note  that,  in  all  cases.  Si  <  Sr-  Let 

us  hrst  treat  the  case  where  Si  >  0,  which  corresponds  to  a  supersonic  flow  from  left 
to  right.  The  Si  wave  would  he  on  the  right  of  the  interface,  and  the  flux  at  a;  =  0, 
F*,  is  then  obviously  given  by  F{Ui).  Similarly,  if  Sr  <  0,  the  flow  is  supersonic,  from 
right  to  left,  and  F*  is  given  by  F{Ur). 

Now,  let  us  examine  the  case  where  Si  <  0  <  Sr-  The  lengths  Xi,Xr  can  then  be 
expressed  as  Xi  =  —TSi  and  Xr  =  TSr-  Expressing  the  integral  form  of  the  Euler 
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equations  on  the  system  in  Fig.  A.l  gives: 

U{x,  0)dx  -  F{U{Xi,  t))dt  +  U{x,  T)dx 
-  j^F{U{Xr,t))dt  +  jl^U{x,Q)dx  =  Q 


(A.3) 


Assuming  piecewise  constant  variables,  and  hence,  piecewise  constant  fluxes,  the  pre¬ 
vious  relation  can  be  re-written  as: 


Ui.{SiT)  -Fi.[T) 

+U\{{Sr  -  Si)T)  -Fr.{-T)+  Ur.{-  A,T)  =  0 


(A.4) 


relation  that,  after  re-arrangement,  leads  to: 

_Fi-  SiUi  -  {Fr  -  SrUr) 

Sr  -Si 


(A.5) 


This  shows  that  once  {Ui,Ur)  is  known  from  the  reconstruction  procedure,  and  once 
the  wavespeeds  (S';,  Sr)  are  estimated,  the  variables  in  the  T*r-region  are  fully  defined. 

The  integral  relation  applied  across  a  given  A;— wave,  k  =  {l,r),  results  in  the 
Rankine-Hugoniot  relations  that  read: 


F*  =  Fi  +  Si  {U*  -  Ui) 

F*  =  Fr  +  Sr  {U*  -  Ur) 


(A.6) 


From  these  2  relations,  one  can  eliminate  U*  in  order  to  determine  F*  as: 

SrFi  -  SiFr  +  SiSriUr  -  U)  ,  . 

Dr  —  Dl 

Thus,  the  flux  evaluated  at  the  i  -|-  1/2  interface  from  the  2-waves  HLL  Riemann 
solver  reads: 

(  Fi  tf  0<Si 

={  F*  tf  S,  <0<  S,  (A.8) 

[  Fr  if  Sr<0 

Once  an  evaluation  of  the  wavespeeds  {Si,  Sr)  is  provided,  the  HLL  Riemann 
solver  is  fully  defined.  A  few  wavespeeds  estimates  can  be  found  in  the  literature 
([13,  18,  19,  52]),  leading  to  schemes  of  different  robustness  and  dissipation.  In  the 
present  development,  the  wavespeeds  are  estimated  following  the  work  of  Einfeldt 
(the  HLL  solver  with  Einfeldt  wave-speeds  is  referred  to  as  HLLE),  as: 

Si  =  min  [qi  —  ci,q  —  c]  Sr  =  max  +  c^,  g  -|-  c]  (A. 9) 

where  U  here  refers  to  the  Roe-averaged  variable  U,  and  c  is  the  local  speed  of  sound. 

This  solver  has  proven  to  be  robust  and  accurate  for  hypersonic  calculations  and 
shock  capturing  purposes.  Its  assumption  of  double  waves  is  however  limiting,  and  the 
consequent  numerical  smearing  of  contact  discontinuities,  shear  waves,  etc...  makes  it 
unsuitable  for  viscous,  turbulent  calculations.  An  extension  of  this  scheme  has  been 
developed  and  presented  by  [53],  where  the  middle  wave  in  the  Riemann  problem  is 
taken  into  account  in  the  derivation  of  the  fluxes.  This  extended  Riemann  solver  is 
named  HLLC  {C  standing  for  Contact),  and  its  derivation  is  given  in  the  next  section. 
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Figure  A. 2:  Three  characteristic  waves,  and  their  characteristic  wavespeeds. 


A. 2  Restoration  of  the  middle  wave  -  the  HLLC 
Riemann  Solver 

This  approximate  Riemann  solver  starts  with  the  assumption  that,  at  a  given  interface 
separating  two  states  in  a  shock  —  tube-like  conhguration,  the  time  evolution  of  a 
system  governed  by  the  Euler  equations  gives  rise  to  three  waves,  of  speed  Si  for  the 
left  moving  wave,  Sr  for  the  right  moving  wave,  and  S*  for  the  contact  wave.  It  is 
then  assumed  that  these  waves  can  be  represented  as  inhnitely  thin  discontinuities, 
separating  constant  states  of  the  fluid.  Si  separates  Ui  from  U^,  S*  is  the  interface 
between  U*  and  U*,  while  Sr  separates  U*  from  Ur-  This  assumption  of  thin  interfaces 
is  justihed  for  both  shocks  and  contact  discontinuities,  and  is  an  approximation  in 
the  case  of  rarefactions.  A  typical  representation  of  a  subsonic  system  is  shown  in 
Fig.  A. 2. 

Let  T  be  the  local  time  step,  T  =  >  0.  Note  that,  in  all  cases.  Si  < 

S*  <  Sr-  Similarly  to  the  HLL  solver,  the  case  of  a  supersonic  flow  from  left  to  right, 
where  Si  >  0,  leads  to  a  flux  at  a;  =  0  given  by  F{Ui).  Again,  if  Sr  <  0,  the  flux  is 
given  by  F{Ur)-  Now,  let  us  examine  the  case  when  Si  <  0  <  Sr-  The  lengths  Xi  and 
Xr  can  then  be  expressed  as  Xi  =  —TSi  and  Xr  =  T Sr  respectively,  and  similarly, 
X*  =  T{S*  —  Si).,  X*  =  T{Sr  —  S*).  The  Euler  equations  in  integral  form,  Eqn.  A.l, 
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can  be  applied  to  the  system  represented  in  Fig.  A. 2,  leading  to  the  relation: 


U(x,0)dx  -  F(U(X,,t))dt  +  U(x,T)dx 
+  fxl'  U(x,T)dx  -  f°F{U(X^,t)}dt  +  f^  U(x,0)dx  ^  0 

With  the  same  assumptions  of  piecewise  constant  variables  and  piecewise  constant 
fluxes  as  in  the  derivation  of  the  HLL  solver,  the  previous  relation  can  be  re-written 
as: 

Ui.{SiT)  -  Fi.{T)  +  -  Si)T) 

+  f/;.  -  S^)T)  -Fr.{-T)+  Ur.  (  -  SrT)  =  0 

relation  that,  after  re-arrangement,  leads  to: 

-  Si)Ut  +  {Sr  -  S*)u:  =  Fi-  SiUi  -  {Fr  -  SrUr) 


(A.ll) 

(A.12) 


This  relation  relates  the  left  and  right  ^-variables  together,  and  is  often  called  the 
consistency  condition.  Note  that  assuming  a  unique  state,  so  that  U*  =  U*  recovers 
the  Consistency  Condition  of  the  HLL  solver,  Eqn.  A. 5. 

The  integral  applied  around  a  control  volume  surrounding  a  given  A;— wave  leads 
to: 

Uk.{SkT)  -  Ffc.(T)  -  Ul.{  -  SkT)  +  Fl{T)  =  0  (A.13) 

so  that  the  corresponding  Rankine-Hugoniot  relations  across  the  /c— wave,  k  =  {l,r), 
are  recovered,  and  reads: 

Fr  =  F,  +  S,  (C/,*  -  O) 

F;  =  Fr  +  Sr  (U;  -  Ur) 

Identically,  the  Rankine-Hugoniot  relation  across  the  T*r-wave  reads: 


f;  =  f;  +  A*  (h;  -  h;)  (A.  15) 

The  relations  A.12,  A. 14  and  A. 15  give  4  relations  for  4  unknowns  (F)*,  F*,  U* ,  and 
U*).  It  is  however  straightforward  to  show  that  they  are  not  linearly  independent.  An 
assumption  has  to  be  made  on  the  intermediate  states  in  order  to  solve  this  system 
of  equations. 

Toro  et  al  [53]  closed  the  relation  by  assuming  that  the  intermediate  wave  had 
the  same  properties  as  a  contact  discontinuity.  Its  propagation  speed  is  then  assumed 
identical  to  the  convective  velocity  in  the  region,  and  this  wave  retains  the  ini¬ 
tial  discontinuity  in  the  passive  scalar  held.  This  also  implies  that  both  convective 
velocities  and  pressures  have  to  match  across  the  interface.  Mathematically,  those 
assumption  are  expressed  by: 

(i?*  .  n  =  g*)  =  {y*  .  n  =  q*)  =  S* 

Pi*  =  P*  (A.  16) 

=  il  ,  0*  =  4>r 

where  0  is  any  passive  scalar  advected  by  the  huid  (0  =  k^^^,Yk,  •••)•  Note  that  one 
special  case  of  passive  scalar  in  the  present  approach  is  the  component  of  the  velocity 
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transverse  to  the  interface,  <pk  =  Vk—  {Vk  ■  nj  n  =  Vk  —  S*n,  k  =  (/,  r).  As  mentioned 

in  Toro  [52],  all  of  these  conditions  are  exactly  satished  by  a  contact  discontinnity 
compnted  from  an  exact  Riemann  solver. 

With  these  assnmptions,  one  can  re-write  the  fonr  first  elements  of  the  vectorial 
eqnation  expressed  in  Eqn.  A.  12,  as: 


{S* 


ur 


u* 


Ur 


Ui 


Si 


Pi 

pI 

Pr 

pi 

P>1 

-  S*) 

pIK 

=  Sr 

PrUr 

-Si 

PlUl 

Plv*l 

p>l 

PrVr 

PlVl 

1 

p*w* 

PrWr 

_  pm  _ 

+ 


pm 

1 - 

piqiui  +  Pin^ 

~h  PrF^x 

piqivi  +  PiUy 

H”  PrFyy 

_  Piqm  +  Pin^  _ 

PrQrWr  +  Pr^tz 

Fi  Fr 


(A.17) 


By  assnmption,  =  S*  for  both  k  =  l/r.  Projecting  the  vectorial  momentnm 
eqnation  on  the  directional  nnit  vector  gives,  along  with  the  hrst  relation,  the  following 
set  of  two  eqnations: 


P*{S*  -  Si)  -h  P*{Sr  -  S*)  =  Pr{Sr  -  Qr)  “  Pl{Sl  -  qi) 
ptiS*  -  Si)  +  ptiSr  -  ^*)1  S*  =  Pl-Pr  +  PrQriSr  -  -  P^Si  -  qi)  (^-18) 


Replacing  the  nnder-braced  term  in  the  last  eqnation  by  the  right-hand  side  of  the 
hrst  eqnation  above  leads  to: 


= 


Pr-Pl+  PiqijSi  -  qi)  -  PrqrjSr  -  qr) 
Pl{Sl  -  qi)  -  Pr{Sr  -  qr) 


(A.19) 


It  is  then  possible  to  manipnlate  Eqn.  A.  14  and  get  the  expressions  for  all  variables. 
One  can  write  the  hrst  fonr  relations,  for  continnity  and  momentnm: 
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Pis* 

PlS*ul  +  P^n, 
PlS*vl  +  P^Uy 
plS*wl  +  P^n, 


PkQk 

PkQk'^k  T  Pk'^x 
PkQk^k  T  Pk^y 
PkQkWk  +  PkUz 
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Pk^k 
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(A.20) 


J 


the  hrst  relation  leads  directly  to  an  expression  for  the  density  in  the  star  region: 

Sk  —  Qk 


Pk  -  Pk 


Sk-S* 


(A.21) 
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Again,  multiplying  the  second  relation  by  n^,  the  third  by  Uy  and  the  last  by  Uz-, 
adding  those  three  relations,  and  using  the  expression  for  p^.  given  in  Eqn.  A. 21  leads 
directly  to: 

PI  =  Pu  +  PkiQk  -  Sk){qk  -  S*)  (A.22) 

Note  that  the  relation  expressed  in  Eqn.  A.22  is  valid  for  both  k  =  (/,  r),  and  satishes 
P*  =  P*.  Let  us  dehne  jSk,  ctk  and  ook  as: 


cifc  —  /dfc  +  1  (A. 23) 

=  —Pk{qk  —  Sk) 

The  state  vectors  can  then  be  expressed  as: 
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—  CnUk  + 


Pk^k^x 

Pk^k^y 

Pk^k^z 

P*S*-Pkqk 

(Sk-S*) 

0 

0 


(A.24) 


With  these  relations,  the  evaluation  of  the  wavespeeds  provides  the  description  of  all 
states  in  the  Riemann  problem  under  consideration.  It  is  then  possible  to  describe 
the  fluxes  at  i  +  1/2  interface.  The  general  expression  for  these  fluxes  is  given  by: 


t^HLLC 

^i+l/2 


Fi 

if 

0<Si 

Fr 

=  Fi  +  Si{U^  -  Ui) 

if 

Si  <0<S* 

f: 

=  Fr  +  Sr{U:  -  Ur) 

if 

S*  <0<  Sr 

Fr 

if 

Sr  <0 

(A.25) 


The  HLLC  Riemann  solver  is  then  fully  dehned,  and  only  the  wavespeeds  are 
needed  to  close  its  expression.  The  estimates  used  in  the  HLLE  solver  are  used  for 
the  evaluation  of  Si  and  whereas  S*  is  dehned  through  relation  A.  19. 
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